Added CustomSourceTime

CHANGELOG.md

@@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Configuration option `config.log_suppression` can be used to control the suppression of log messages.
 - `abort()` for `Job` and `mode solver`, Job or mode solver whose status is not success or error(e.g. running, draft) can be aborted, if Job or mode solver is abort, it can't be submitted, a new one needs to be created and submitted.
 - `web.abort()` and `Job.abort()` methods allowing to abort running tasks without deleting them. If a task is aborted, it cannot be restarted later, a new one needs to be created and submitted.
+- Source with arbitrary user-specified time dependence through `CustomSourceTime`.
 
 ### Changed
 - Add width and height options to Simulation.plot_3d

tests/sims/simulation_2_3_0.json

@@ -929,6 +929,34 @@
             "angle_phi": 0.6283185307179586,
             "pol_angle": 0.0,
             "injection_axis": 2
+        },
+        {
+            "type": "UniformCurrentSource",
+            "center": [
+                0.0,
+                0.5,
+                0.0
+            ],
+            "size": [
+                0.0,
+                0.0,
+                0.0
+            ],
+            "source_time": {
+                "amplitude": 1.0,
+                "phase": 0.0,
+                "type": "CustomSourceTime",
+                "freq0": 200000000000000.0,
+                "fwidth": 40000000000000.0,
+                "offset": 5.0,
+                "source_time_dataset": {
+                    "type": "TimeDataset",
+                    "values": "TimeDataArray"
+                }
+            },
+            "name": null,
+            "interpolate": true,
+            "polarization": "Hx"
         }
     ],
     "boundary_spec": {

tests/sims/simulation_2_4_0rc1.json

@@ -960,6 +960,34 @@
             "angle_phi": 0.6283185307179586,
             "pol_angle": 0.0,
             "injection_axis": 2
+        },
+        {
+            "type": "UniformCurrentSource",
+            "center": [
+                0.0,
+                0.5,
+                0.0
+            ],
+            "size": [
+                0.0,
+                0.0,
+                0.0
+            ],
+            "source_time": {
+                "amplitude": 1.0,
+                "phase": 0.0,
+                "type": "CustomSourceTime",
+                "freq0": 200000000000000.0,
+                "fwidth": 40000000000000.0,
+                "offset": 5.0,
+                "source_time_dataset": {
+                    "type": "TimeDataset",
+                    "values": "TimeDataArray"
+                }
+            },
+            "name": null,
+            "interpolate": true,
+            "polarization": "Hx"
         }
     ],
     "boundary_spec": {

tests/test_components/test_IO.py

@@ -37,6 +37,8 @@ def set_datasets_to_none(sim):
             src["field_dataset"] = None
         elif src["type"] == "CustomCurrentSource":
             src["current_dataset"] = None
+        if src["source_time"]["type"] == "CustomSourceTime":
+            src["source_time"]["source_time_dataset"] = None
     for structure in sim_dict["structures"]:
         if structure["geometry"]["type"] == "TriangleMesh":
             structure["geometry"]["mesh_dataset"] = None

tests/test_components/test_source.py

@@ -13,6 +13,9 @@
 S = td.PointDipole(source_time=ST, polarization="Ex")
 
 
+ATOL = 1e-8
+
+
 def test_plot_source_time():
 
     for val in ("real", "imag", "abs"):
@@ -247,3 +250,41 @@ def check_freq_grid(freq_grid, num_freqs):
             mode_index=0,
             num_freqs=-10,
         )
+
+
+def test_custom_source_time():
+    g = td.GaussianPulse(freq0=1, fwidth=0.1)
+    ts = np.linspace(0, 30, 1001)
+    amp_time = g.amp_time(ts)
+
+    # basic test
+    cst = td.CustomSourceTime.from_values(freq0=1, fwidth=0.1, values=amp_time, dt=ts[1] - ts[0])
+    assert np.allclose(cst.amp_time(ts), amp_time, rtol=0, atol=ATOL)
+
+    # test single value validation error
+    with pytest.raises(pydantic.ValidationError):
+        vals = td.components.data.data_array.TimeDataArray([1], coords=dict(t=[0]))
+        dataset = td.components.data.dataset.TimeDataset(values=vals)
+        cst = td.CustomSourceTime(source_time_dataset=dataset, freq0=1, fwidth=0.1)
+        assert np.allclose(cst.amp_time([0]), [1], rtol=0, atol=ATOL)
+
+    # test interpolation
+    cst = td.CustomSourceTime.from_values(freq0=1, fwidth=0.1, values=np.linspace(0, 9, 10), dt=0.1)
+    assert np.allclose(cst.amp_time(0.09), [0.9], rtol=0, atol=ATOL)
+
+    # test sampling warning
+    cst = td.CustomSourceTime.from_values(freq0=1, fwidth=0.1, values=np.linspace(0, 9, 10), dt=0.1)
+    source = td.PointDipole(center=(0, 0, 0), source_time=cst, polarization="Ex")
+    sim = td.Simulation(
+        size=(10, 10, 10),
+        run_time=1e-12,
+        grid_spec=td.GridSpec.uniform(dl=0.1),
+        sources=[source],
+    )
+
+    # test out of range validation error
+    with pytest.raises(td.exceptions.ValidationError):
+        vals = np.cos(sim.tmesh[:-1])
+        cst = td.CustomSourceTime.from_values(freq0=1, fwidth=0.1, values=vals, dt=sim.dt)
+        source = td.PointDipole(center=(0, 0, 0), source_time=cst, polarization="Ex")
+        sim = sim.updated_copy(sources=[source])

tests/utils.py

@@ -365,6 +365,14 @@ def prepend_tmp(path):
             angle_phi=np.pi / 5,
             injection_axis=2,
         ),
+        UniformCurrentSource(
+            size=(0, 0, 0),
+            center=(0, 0.5, 0),
+            polarization="Hx",
+            source_time=CustomSourceTime.from_values(
+                freq0=2e14, fwidth=4e13, values=np.linspace(0, 10, 1000), dt=1e-12 / 100
+            ),
+        ),
     ],
     monitors=(
         FieldMonitor(

tidy3d/__init__.py

@@ -26,7 +26,7 @@
 from .components.apodization import ApodizationSpec
 
 # sources
-from .components.source import GaussianPulse, ContinuousWave
+from .components.source import GaussianPulse, ContinuousWave, CustomSourceTime
 from .components.source import UniformCurrentSource, PlaneWave, ModeSource, PointDipole
 from .components.source import GaussianBeam, AstigmaticGaussianBeam
 from .components.source import CustomFieldSource, TFSF, CustomCurrentSource

tidy3d/components/data/dataset.py

@@ -12,6 +12,7 @@
 from .data_array import ScalarFieldDataArray, ScalarFieldTimeDataArray, ScalarModeFieldDataArray
 from .data_array import ModeIndexDataArray
 from .data_array import TriangleMeshDataArray
+from .data_array import TimeDataArray
 
 from ..base import Tidy3dBaseModel
 from ..types import Axis
@@ -403,3 +404,11 @@ class TriangleMeshDataset(Dataset):
         description="Dataset containing the surface triangles and corresponding face indices "
         "for a surface mesh.",
     )
+
+
+class TimeDataset(Dataset):
+    """Dataset for storing a function of time."""
+
+    values: TimeDataArray = pd.Field(
+        ..., title="Values", description="Values as a function of time."
+    )

tidy3d/components/simulation.py

@@ -25,7 +25,7 @@
 from .boundary import PML, StablePML, Absorber, AbsorberSpec
 from .structure import Structure
 from .source import SourceType, PlaneWave, GaussianBeam, AstigmaticGaussianBeam, CustomFieldSource
-from .source import CustomCurrentSource
+from .source import CustomCurrentSource, CustomSourceTime
 from .source import TFSF, Source
 from .monitor import MonitorType, Monitor, FreqMonitor, SurfaceIntegrationMonitor
 from .monitor import AbstractFieldMonitor, DiffractionMonitor, AbstractFieldProjectionMonitor
@@ -839,6 +839,7 @@ def _post_init_validators(self) -> None:
         """Call validators taking z`self` that get run after init."""
         self._validate_no_structures_pml()
         self._validate_tfsf_nonuniform_grid()
+        self._validate_customsourcetime()
 
     def _validate_no_structures_pml(self) -> None:
         """Ensure no structures terminate / have bounds inside of PML."""
@@ -908,6 +909,31 @@ def _validate_tfsf_nonuniform_grid(self) -> None:
                             f"axis, '{'xyz'[source.injection_axis]}'."
                         )
 
+    def _validate_customsourcetime(self) -> None:
+        """Make sure custom source time is not undersampled.
+        Also, make sure that all simulation.tmesh values are covered."""
+        for source in self.sources:
+            if isinstance(source.source_time, CustomSourceTime):
+                dataset = source.source_time.source_time_dataset
+                if dataset is None:
+                    continue
+                times = dataset.values.coords["t"].values
+                if min(times) > self.tmesh[0] or max(times) < self.tmesh[-1]:
+                    raise ValidationError(
+                        "'CustomSourceTime' found with time coordinates "
+                        "'times' that do not cover the entire 'Simulation.tmesh'. Currently, "
+                        f"'(min(times), max(times)) = ({min(times)}, {max(times)})', while "
+                        f"'(min(tmesh), max(tmesh)) = ({self.tmesh[0]}, {self.tmesh[-1]}).' "
+                    )
+                max_dt = np.amax(np.diff(times))
+                if max_dt > self.dt * 1.01:
+                    log.warning(
+                        f"'CustomSourceTime' found with time step 'max(dt) = {max_dt:.3g}', "
+                        f"while the simulation time step is 'dt={self.dt}'. "
+                        "We recommend that the largest time step of the custom source "
+                        f"be smaller than the time step of the simulation."
+                    )
+
     """ Pre submit validation (before web.upload()) """
 
     def validate_pre_upload(self, source_required: bool = True) -> None:
@@ -2683,6 +2709,11 @@ def custom_datasets(self) -> List[Dataset]:
         """List of custom datasets for verification purposes. If the list is not empty, then
         the simulation needs to be exported to hdf5 to store the data.
         """
+        datasets_source_time = [
+            src.source_time.source_time_dataset
+            for src in self.sources
+            if isinstance(src.source_time, CustomSourceTime)
+        ]
         datasets_field_source = [
             src.field_dataset for src in self.sources if isinstance(src, CustomFieldSource)
         ]
@@ -2699,7 +2730,13 @@ def custom_datasets(self) -> List[Dataset]:
                 for geometry in struct.geometry.geometries:
                     datasets_geometry += geometry.mesh_dataset
 
-        return datasets_field_source + datasets_current_source + datasets_medium + datasets_geometry
+        return (
+            datasets_source_time
+            + datasets_field_source
+            + datasets_current_source
+            + datasets_medium
+            + datasets_geometry
+        )
 
     def _volumetric_structures_grid(self, grid: Grid) -> Tuple[Structure]:
         """Generate a tuple of structures wherein any 2D materials are converted to 3D

tidy3d/components/source.py

@@ -1,5 +1,7 @@
 """Defines electric current sources for injecting light into simulation."""
+# pylint: disable=too-many-lines
 
+from __future__ import annotations
 from abc import ABC, abstractmethod
 from typing import Union, Tuple, Optional
 
@@ -9,17 +11,19 @@
 
 from .base import Tidy3dBaseModel, cached_property
 
-from .types import Direction, Polarization, Ax, FreqBound, ArrayFloat1D, Axis, PlotVal
+from .types import Direction, Polarization, Ax, FreqBound
+from .types import ArrayFloat1D, Axis, PlotVal, ArrayComplex1D
 from .validators import assert_plane, assert_volumetric, validate_name_str, get_value
 from .validators import warn_if_dataset_none, assert_single_freq_in_range
-from .data.dataset import FieldDataset
+from .data.dataset import FieldDataset, TimeDataset
+from .data.data_array import TimeDataArray
 from .geometry import Box, Coordinate
 from .mode import ModeSpec
 from .viz import add_ax_if_none, PlotParams, plot_params_source
 from .viz import ARROW_COLOR_SOURCE, ARROW_ALPHA, ARROW_COLOR_POLARIZATION
 from ..constants import RADIAN, HERTZ, MICROMETER, GLANCING_CUTOFF
 from ..constants import inf  # pylint:disable=unused-import
-from ..exceptions import SetupError
+from ..exceptions import SetupError, ValidationError
 from ..log import log
 
 
@@ -291,7 +295,90 @@ def amp_time(self, time: float) -> complex:
         return const * offset * oscillation * amp
 
 
-SourceTimeType = Union[GaussianPulse, ContinuousWave]
+class CustomSourceTime(Pulse):
+    """Custom source time dependence, real or complex valued.
+
+    Note
+    ----
+    The source time dependence is linearly interpolated to the simulation time steps.
+    To ensure that this interpolation does not introduce artifacts, it is necessary
+    to use a sampling rate that is sufficiently fast relative to the simulation time step.
+
+    Example
+    -------
+    >>> cst = td.CustomSourceTime.from_values(freq0=1, fwidth=0.1,
+    >>>     values=np.linspace(0, 9, 10), dt=0.1)
+    """
+
+    source_time_dataset: Optional[TimeDataset] = pydantic.Field(
+        ..., title="Source time dataset", description="Dataset for storing the custom source time."
+    )
+
+    _source_time_dataset_none_warning = warn_if_dataset_none("source_time_dataset")
+
+    @pydantic.validator("source_time_dataset", always=True)
+    def _more_than_one_time(cls, val):
+        """Must have more than one time to interpolate."""
+        if val is None:
+            return val
+        if val.values.size <= 1:
+            raise ValidationError("'CustomSourceTime' must have more than one time coordinate.")
+        return val
+
+    @classmethod
+    def from_values(
+        cls, freq0: float, fwidth: float, values: ArrayComplex1D, dt: float
+    ) -> CustomSourceTime:
+        """Create a :class:`.CustomSourceTime` from a numpy array.
+
+        Parameters
+        ----------
+        freq0 : float
+            Estimated central frequency of the source.
+        fwidth : float
+            Estimated frequency width of the source.
+        values: ArrayComplex1D
+            Complex values of the source amplitude.
+        dt: float
+            Time step for the `values` array. This value should be sufficiently small
+            relative to the simulation `dt` in order to avoid interpolation artifacts.
+
+
+        Returns
+        -------
+        CustomSourceTime
+            :class:`.CustomSourceTime` with these values, and time coordinates evenly spaced
+            between 0 and dt * (N-1) with a step size of `dt`, where N is the length of
+            the values array.
+        """
+        times = np.arange(len(values)) * dt
+        source_time_dataarray = TimeDataArray(values, coords=dict(t=times))
+        source_time_dataset = TimeDataset(values=source_time_dataarray)
+        return CustomSourceTime(
+            freq0=freq0,
+            fwidth=fwidth,
+            source_time_dataset=source_time_dataset,
+        )
+
+    def amp_time(self, time: float) -> complex:
+        """Complex-valued source amplitude as a function of time.
+
+        Parameters
+        ----------
+        time : float
+            Time in seconds.
+
+        Returns
+        -------
+        complex
+            Complex-valued source amplitude at that time.
+        """
+        if self.source_time_dataset is None:
+            return None
+        return self.source_time_dataset.values.interp(t=time).to_numpy()
+
+
+SourceTimeType = Union[GaussianPulse, ContinuousWave, CustomSourceTime]
 
 """ Source objects """