Revert to num_freqs=1 for broadband angled gaussian beam

momchil-flex · momchil-flex · commit c70cba7e6042 · 2025-06-25T21:45:10.000+02:00
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -18,6 +18,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Arrow lengths are now scaled consistently in the X and Y directions, and their lengths no longer exceed the height of the plot window.
 - Bug in `PlaneWave` defined with a negative `angle_theta` which would lead to wrong injection.
 
+### Changed
+- `GaussianBeam` and `AstigmaticGaussianBeam` default `num_freqs` reset to 1 (it was set to 3 in v2.8.0) and a warning is issued for a broadband, angled beam for which `num_freqs` may not be sufficiently large.
+- Set the maximum `num_freqs` to 20 for all broadband sources (we have been warning about the introduction of this hard limit for a while).
+
 ## [2.9.0rc1] - 2025-06-10
 
 ### Added
diff --git a/tests/test_components/test_source.py b/tests/test_components/test_source.py
@@ -459,3 +459,80 @@ def test_fixed_angle_source():
     )
 
     assert not plane_wave._is_fixed_angle
+
+
+def test_broadband_angled_gaussian_warning():
+    g = td.GaussianPulse(freq0=1e14, fwidth=0.8e14)
+    # Case 1: num_freqs = 3, angle_theta = np.pi / 3, should warn
+    with AssertLogLevel("WARNING", contains_str="number of frequencies"):
+        s = td.GaussianBeam(
+            size=(0, 1, 1),
+            source_time=g,
+            pol_angle=np.pi / 2,
+            direction="+",
+            angle_theta=np.pi / 3,
+            num_freqs=3,
+        )
+        _ = td.Simulation(
+            size=(2, 2, 2),
+            run_time=1e-12,
+            grid_spec=td.GridSpec.uniform(dl=0.1),
+            sources=[s],
+            normalize_index=None,
+        )
+
+    # Case 2: Increasing to num_freqs = 10 should NOT warn
+    with AssertLogLevel(None):
+        s = td.GaussianBeam(
+            size=(0, 1, 1),
+            source_time=g,
+            pol_angle=np.pi / 2,
+            direction="+",
+            angle_theta=np.pi / 3,
+            num_freqs=10,
+        )
+        _ = td.Simulation(
+            size=(2, 2, 2),
+            run_time=1e-12,
+            grid_spec=td.GridSpec.uniform(dl=0.1),
+            sources=[s],
+            normalize_index=None,
+        )
+
+    # Case 3: Case 2 but changed to astigmatic gaussian beam with one larger waist size should warn
+    with AssertLogLevel("WARNING", contains_str="number of frequencies"):
+        s = td.AstigmaticGaussianBeam(
+            size=(0, 1, 1),
+            source_time=g,
+            pol_angle=np.pi / 2,
+            direction="+",
+            angle_theta=np.pi / 3,
+            num_freqs=10,
+            waist_sizes=(1, 5),
+        )
+        _ = td.Simulation(
+            size=(2, 2, 2),
+            run_time=1e-12,
+            grid_spec=td.GridSpec.uniform(dl=0.1),
+            sources=[s],
+            normalize_index=None,
+        )
+
+    # Case 4: Case 3 but with num_freqs = 1 should NOT warn (broadband treatment is off)
+    with AssertLogLevel(None):
+        s = td.AstigmaticGaussianBeam(
+            size=(0, 1, 1),
+            source_time=g,
+            pol_angle=np.pi / 2,
+            direction="+",
+            angle_theta=np.pi / 3,
+            num_freqs=1,
+            waist_sizes=(1, 5),
+        )
+        _ = td.Simulation(
+            size=(2, 2, 2),
+            run_time=1e-12,
+            grid_spec=td.GridSpec.uniform(dl=0.1),
+            sources=[s],
+            normalize_index=None,
+        )
diff --git a/tests/test_package/test_log.py b/tests/test_package/test_log.py
@@ -85,13 +85,13 @@ def test_logging_warning_capture():
         name="mode",
     )
 
-    # 2 warnings: too high num_freqs; too many points
+    # 1 warning: too many points
     mode_source = td.ModeSource(
         size=(domain_size, 0, domain_size),
         source_time=source_time,
         mode_spec=td.ModeSpec(num_modes=2, precision="single"),
         mode_index=1,
-        num_freqs=50,
+        num_freqs=10,
         direction="-",
     )
 
@@ -218,7 +218,7 @@ def test_logging_warning_capture():
     sim.validate_pre_upload()
     warning_list = td.log.captured_warnings()
     print(json.dumps(warning_list, indent=4))
-    assert len(warning_list) == 30
+    assert len(warning_list) == 29
     td.log.set_capture(False)
 
     # check that capture doesn't change validation errors
diff --git a/tidy3d/components/simulation.py b/tidy3d/components/simulation.py
@@ -3618,6 +3618,53 @@ def _source_homogeneous_isotropic(cls, val, values):
                                 "A fixed angle plane wave can only be injected into a homogeneous isotropic"
                                 "dispersionless medium."
                             )
+                    # check if broadband angled gaussian beam frequency variation is too fast
+                    if (
+                        isinstance(source, (GaussianBeam, AstigmaticGaussianBeam))
+                        and np.abs(source.angle_theta) > 0
+                        and source.num_freqs > 1
+                    ):
+
+                        def radius(waist_radius, waist_distance, k0):
+                            """Gaussian beam radius at a given waist distance and k0."""
+                            z_r = waist_radius**2 * k0 / 2
+                            return waist_radius * np.sqrt(1 + (waist_distance / z_r) ** 2)
+
+                        # A slanted GaussianBeam will accumulate a phase that's frequency-dependent
+                        # like phi = K f, with the derivative dphi / df = K = 2 * pi * n * r * sin(theta) / c_0.
+                        # Here, we compute the maximum value of this coefficient computed at the waist radius
+                        # and over all frequencies. Then we compare this to the frequency spacing to
+                        # determine whether the frequency dependence is too fast, and issue a warning.
+                        optical_path_length = []
+                        freqs = source.frequency_grid
+                        for freq in freqs:
+                            n_freq, _ = src_medium.nk_model(frequency=freq)
+                            k0 = 2 * np.pi * n_freq * freq / C_0
+                            if isinstance(source, GaussianBeam):
+                                rad = radius(source.waist_radius, source.waist_distance, k0)
+                            else:
+                                rad = max(
+                                    radius(source.waist_sizes[0], source.waist_distances[0], k0),
+                                    radius(source.waist_sizes[1], source.waist_distances[1], k0),
+                                )
+                            optical_path_length.append(n_freq * rad * np.sin(source.angle_theta))
+                        # Maximum value of the path length over all freqs
+                        max_path_length = np.max(optical_path_length)
+                        # Maximum value of the phase difference
+                        max_phase_diff = max_path_length * 2 * np.pi * (freqs[-1] - freqs[0]) / C_0
+                        # Compare this in magnitude to the frequency spacing assuming uniform
+                        # spacing. This is heuristic since in reality we use a Chebyshev grid,
+                        # but it should be a good rule of thumb. Because the Chebyshev interpolation
+                        # is much better than simple interpolation, we don't require << 1, just < 1
+                        if not max_phase_diff / source.num_freqs < 1:
+                            log.warning(
+                                f"Broadband, angled {source.type} source has a phase dependence "
+                                "with frequency that might be under-resolved by the provided "
+                                "number of frequencies. Consider reducing the source bandwidth, "
+                                "or increasing the 'num_freqs' of the source, and verify the "
+                                "source injection in an empty simulation.",
+                            )
+
         return val
 
     @pydantic.validator("normalize_index", always=True)
diff --git a/tidy3d/components/source/field.py b/tidy3d/components/source/field.py
@@ -27,8 +27,6 @@
 
 # width of Chebyshev grid used for broadband sources (in units of pulse width)
 CHEB_GRID_WIDTH = 1.5
-# Number of frequencies in a broadband source above which to issue a warning
-WARN_NUM_FREQS = 20
 # For broadband plane waves with constan in-plane k, the Chebyshev grid is truncated at
 # ``CRITICAL_FREQUENCY_FACTOR * f_crit``, where ``f_crit`` is the critical frequency
 # (oblique propagation).
@@ -90,14 +88,14 @@ def _dir_vector(self) -> tuple[float, float, float]:
 class BroadbandSource(Source, ABC):
     """A source with frequency dependent field distributions."""
 
-    # Default as for analytic beam sources; overwrriten for ModeSource below
     num_freqs: int = pydantic.Field(
-        3,
+        1,
         title="Number of Frequency Points",
         description="Number of points used to approximate the frequency dependence of the injected "
-        "field. Default is 3, which should cover even very broadband sources. For simulations "
-        "which are not very broadband and the source is very large (e.g. metalens simulations), "
-        "decreasing the value to 1 may lead to a speed up in the preprocessing.",
+        "field. A Chebyshev interpolation is used, thus, only a small number of points is "
+        "typically sufficient to obtain converged results. Note that larger values of 'num_freqs' "
+        "could spread out the source time signal and introduce numerical noise, or prevent timely  "
+        "field decay.",
         ge=1,
         le=20,
     )
@@ -116,23 +114,6 @@ def _chebyshev_freq_grid(self, freq_min, freq_max):
         cheb_points = np.cos(np.pi * np.flip(uni_points))
         return freq_avg + freq_diff * cheb_points
 
-    @pydantic.validator("num_freqs", always=True, allow_reuse=True)
-    def _warn_if_large_number_of_freqs(cls, val):
-        """Warn if a large number of frequency points is requested."""
-
-        if val is None:
-            return val
-
-        if val >= WARN_NUM_FREQS:
-            log.warning(
-                f"A large number ({val}) of frequency points is used in a broadband source. "
-                "This can lead to solver slow-down and increased cost, and even introduce "
-                "numerical noise. This may become a hard limit in future Tidy3D versions.",
-                custom_loc=["num_freqs"],
-            )
-
-        return val
-
 
 """ Source current profiles determined by user-supplied data on a plane."""
 
@@ -422,16 +403,6 @@ class ModeSource(DirectionalSource, PlanarSource, BroadbandSource):
         "``num_modes`` in the solver will be set to ``mode_index + 1``.",
     )
 
-    num_freqs: int = pydantic.Field(
-        1,
-        title="Number of Frequency Points",
-        description="Number of points used to approximate the frequency dependence of injected "
-        "field. A Chebyshev interpolation is used, thus, only a small number of points, i.e., less "
-        "than 20, is typically sufficient to obtain converged results.",
-        ge=1,
-        le=99,
-    )
-
     @cached_property
     def angle_theta(self):
         """Polar angle of propagation."""
@@ -515,6 +486,17 @@ class PlaneWave(AngledFieldSource, PlanarSource, BroadbandSource):
         discriminator=TYPE_TAG_STR,
     )
 
+    num_freqs: int = pydantic.Field(
+        3,
+        title="Number of Frequency Points",
+        description="Number of points used to approximate the frequency dependence of the injected "
+        "field. Default is 3, which should cover even very broadband plane waves. For simulations "
+        "which are not very broadband and the source is very large (e.g. metalens simulations), "
+        "decreasing the value to 1 may lead to a speed up in the preprocessing.",
+        ge=1,
+        le=20,
+    )
+
     @cached_property
     def _is_fixed_angle(self) -> bool:
         """Whether the plane wave is at a fixed non-zero angle."""
@@ -594,6 +576,18 @@ class GaussianBeam(AngledFieldSource, PlanarSource, BroadbandSource):
         units=MICROMETER,
     )
 
+    num_freqs: int = pydantic.Field(
+        1,
+        title="Number of Frequency Points",
+        description="Number of points used to approximate the frequency dependence of the injected "
+        "field. For broadband, angled Gaussian beams it is advisable to check the beam propagation "
+        "in an empty simulation to ensure there are no injection artifacts when 'num_freqs' > 1. "
+        "Note that larger values of 'num_freqs' could spread out the source time signal and "
+        "introduce numerical noise, or prevent timely field decay.",
+        ge=1,
+        le=20,
+    )
+
 
 class AstigmaticGaussianBeam(AngledFieldSource, PlanarSource, BroadbandSource):
     """The simple astigmatic Gaussian distribution allows
@@ -642,6 +636,18 @@ class AstigmaticGaussianBeam(AngledFieldSource, PlanarSource, BroadbandSource):
         units=MICROMETER,
     )
 
+    num_freqs: int = pydantic.Field(
+        1,
+        title="Number of Frequency Points",
+        description="Number of points used to approximate the frequency dependence of the injected "
+        "field. For broadband, angled Gaussian beams it is advisable to check the beam propagation "
+        "in an empty simulation to ensure there are no injection artifacts when 'num_freqs' > 1. "
+        "Note that larger values of 'num_freqs' could spread out the source time signal and "
+        "introduce numerical noise, or prevent timely field decay.",
+        ge=1,
+        le=20,
+    )
+
 
 class TFSF(AngledFieldSource, VolumeSource, BroadbandSource):
     """Total-field scattered-field (TFSF) source that can inject a plane wave in a finite region.
