From 00f56d254896717a7625d092d6992acf81bb54ae Mon Sep 17 00:00:00 2001 From: Casey Wojcik Date: Fri, 19 May 2023 11:04:15 +0100 Subject: [PATCH] Added CustomSourceTime --- CHANGELOG.md | 1 + tests/sims/simulation_2_3_0.json | 28 ++++++++ tests/sims/simulation_2_4_0rc1.json | 28 ++++++++ tests/test_components/test_IO.py | 2 + tests/test_components/test_source.py | 41 ++++++++++++ tests/utils.py | 8 +++ tidy3d/__init__.py | 2 +- tidy3d/components/data/dataset.py | 9 +++ tidy3d/components/simulation.py | 44 ++++++++++++- tidy3d/components/source.py | 95 ++++++++++++++++++++++++++-- 10 files changed, 251 insertions(+), 7 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6ff246178e..c57d447374 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/tests/sims/simulation_2_3_0.json b/tests/sims/simulation_2_3_0.json index 6a5a1fb6af..609b27d1a4 100644 --- a/tests/sims/simulation_2_3_0.json +++ b/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": { diff --git a/tests/sims/simulation_2_4_0rc1.json b/tests/sims/simulation_2_4_0rc1.json index e8aa565e0c..fe7610f04e 100644 --- a/tests/sims/simulation_2_4_0rc1.json +++ b/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": { diff --git a/tests/test_components/test_IO.py b/tests/test_components/test_IO.py index f4ccf7e464..66a0fefe42 100644 --- a/tests/test_components/test_IO.py +++ b/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 diff --git a/tests/test_components/test_source.py b/tests/test_components/test_source.py index fd2d18b63c..a053688298 100644 --- a/tests/test_components/test_source.py +++ b/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]) diff --git a/tests/utils.py b/tests/utils.py index b4790bfab4..470c76a5ba 100644 --- a/tests/utils.py +++ b/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( diff --git a/tidy3d/__init__.py b/tidy3d/__init__.py index a221b434db..084354cb0d 100644 --- a/tidy3d/__init__.py +++ b/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 diff --git a/tidy3d/components/data/dataset.py b/tidy3d/components/data/dataset.py index 939df3bd2d..3d80ca4cb4 100644 --- a/tidy3d/components/data/dataset.py +++ b/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." + ) diff --git a/tidy3d/components/simulation.py b/tidy3d/components/simulation.py index 6c55032373..54d6274bab 100644 --- a/tidy3d/components/simulation.py +++ b/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,9 @@ # 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 +CUSTOMSOURCETIME_TOLERANCE = 1.01 + class Simulation(Box): # pylint:disable=too-many-public-methods """Contains all information about Tidy3d simulation. @@ -839,6 +842,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 +912,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 * CUSTOMSOURCETIME_TOLERANCE: + 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 +2712,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 +2733,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 diff --git a/tidy3d/components/source.py b/tidy3d/components/source.py index f8de1836bc..1150b79f3f 100644 --- a/tidy3d/components/source.py +++ b/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 """