Added CustomSourceTime

CHANGELOG.md

@@ -10,6 +10,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - `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.
 - `FastDispersionFitter` for fast fitting of material dispersion data.
+- 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,46 @@ 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[:-2])
+        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])
+
+    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
@@ -68,6 +68,10 @@
 # height of the PML plotting boxes along any dimensions where sim.size[dim] == 0
 PML_HEIGHT_FOR_0_DIMS = 0.02
 
+# allow some numerical flexibility before warning about CustomSourceTime
+# in units of dt
+CUSTOMSOURCETIME_TOL = 1.1
+
 
 class Simulation(Box):  # pylint:disable=too-many-public-methods
     """Contains all information about Tidy3d simulation.
@@ -839,6 +843,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 +913,34 @@ 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] - CUSTOMSOURCETIME_TOL * self.dt
+                ):
+                    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 * CUSTOMSOURCETIME_TOL:
+                    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 +2716,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 +2737,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