Changed flux and DFT definitions with complex fields to only use real part, and fixed TFSF flux leakage bug

Author: caseyflex

CHANGELOG.md

@@ -16,8 +16,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Credit cost for remote mode solver has been modified to be defined in advance based on the mode solver details. Previously, the cost was based on elapsed runtime. On average, there should be little difference in the cost.
 - Mode solves that are part of an FDTD simulation (i.e. for mode sources and monitors) are now charged at the same flex credit cost as a corresponding standalone mode solver call.
 - Any `FreqMonitor.freqs` or `Source.source_time.freq0` smaller than `1e5` now raise an error as this must be incorrect setup that is outside the Tidy3D intended range (note default frequency is `Hz`).
+- When using complex fields (e.g. with Bloch boundaries), FluxTimeMonitor and frequency-domain fields (including derived quantities like flux) now only use the real part of the time-domain electric field.
 
 ### Fixed
+- Fixed energy leakage in TFSF when using complex fields.
 
 ## [2.5.0rc3] - 2023-11-30
 

tests/sims/simulation_2_5_0rc3.h5

@@ -843,8 +843,8 @@ def poynting(self) -> ScalarFieldTimeDataArray:
         # Tangential fields are ordered as E1, E2, H1, H2
         tan_fields = self._colocated_tangential_fields
         dim1, dim2 = self._tangential_dims
-        e_x_h = tan_fields["E" + dim1] * tan_fields["H" + dim2]
-        e_x_h -= tan_fields["E" + dim2] * tan_fields["H" + dim1]
+        e_x_h = np.real(tan_fields["E" + dim1]) * np.real(tan_fields["H" + dim2])
+        e_x_h -= np.real(tan_fields["E" + dim2]) * np.real(tan_fields["H" + dim1])
         return e_x_h
 
     @cached_property

tidy3d/components/data/monitor_data.py

@@ -11,8 +11,6 @@
 
 from .monitor_data import MonitorDataTypes, MonitorDataType, AbstractFieldData, FieldTimeData
 from ..simulation import Simulation
-from ..boundary import BlochBoundary
-from ..source import TFSF
 from ..types import Ax, Axis, annotate_type, FieldVal, PlotScale, ColormapType
 from ..viz import equal_aspect, add_ax_if_none
 from ...exceptions import DataError, Tidy3dKeyError
@@ -116,16 +114,10 @@ def source_spectrum(self, source_index: int) -> Callable:
         times = self.simulation.tmesh
         dt = self.simulation.dt
 
-        # get boundary information to determine whether to use complex fields
-        boundaries = self.simulation.boundary_spec.to_list
-        boundaries_1d = [boundary_1d for dim_boundary in boundaries for boundary_1d in dim_boundary]
-        complex_fields = any(isinstance(boundary, BlochBoundary) for boundary in boundaries_1d)
-        complex_fields = complex_fields and not isinstance(source, TFSF)
-
         # plug in mornitor_data frequency domain information
         def source_spectrum_fn(freqs):
             """Source amplitude as function of frequency."""
-            spectrum = source_time.spectrum(times, freqs, dt, complex_fields)
+            spectrum = source_time.spectrum(times, freqs, dt)
 
             # Remove user defined amplitude and phase from the normalization
             # such that they would still have an effect on the output fields.

tidy3d/components/data/sim_data.py

@@ -155,6 +155,11 @@ def _get_n0(
             )
         return n0
 
+    @property
+    def complex_fields(self) -> bool:
+        """Whether the model uses complex fields."""
+        pass
+
 
 class NonlinearSusceptibility(NonlinearModel):
     """Model for an instantaneous nonlinear chi3 susceptibility.
@@ -215,6 +220,11 @@ def _validate_numiters(cls, val):
             )
         return val
 
+    @property
+    def complex_fields(self) -> bool:
+        """Whether the model uses complex fields."""
+        return False
+
 
 class TwoPhotonAbsorption(NonlinearModel):
     """Model for two-photon absorption (TPA) nonlinearity which gives an intensity-dependent
@@ -287,6 +297,11 @@ def _validate_medium(self, medium: AbstractMedium):
         if self.n0 is not None:
             self._validate_medium_freqs(medium, [])
 
+    @property
+    def complex_fields(self) -> bool:
+        """Whether the model uses complex fields."""
+        return True
+
 
 class KerrNonlinearity(NonlinearModel):
     """Model for Kerr nonlinearity which gives an intensity-dependent refractive index
@@ -358,6 +373,11 @@ def _validate_medium(self, medium: AbstractMedium):
         if self.n0 is not None:
             self._validate_medium_freqs(medium, [])
 
+    @property
+    def complex_fields(self) -> bool:
+        """Whether the model uses complex fields."""
+        return True
+
 
 NonlinearModelType = Union[NonlinearSusceptibility, TwoPhotonAbsorption, KerrNonlinearity]
 

tidy3d/components/medium.py

@@ -2403,6 +2403,10 @@ def complex_fields(self) -> bool:
         """
         if any(isinstance(boundary[0], BlochBoundary) for boundary in self.boundary_spec.to_list):
             return True
+        for medium in self.scene.mediums:
+            if medium.nonlinear_spec is not None:
+                if any(model.complex_fields for model in medium._nonlinear_models):
+                    return True
         return False
 
     @cached_property

tidy3d/components/simulation.py

@@ -46,9 +46,9 @@ def plot_spectrum(
         num_freqs: int = 101,
         val: PlotVal = "real",
         ax: Ax = None,
-        complex_fields: bool = False,
     ) -> Ax:
         """Plot the complex-valued amplitude of the source time-dependence.
+        Note: Only the real part of the time signal is used.
 
         Parameters
         ----------
@@ -61,8 +61,6 @@ def plot_spectrum(
             Number of frequencies to plot within the SourceTime.frequency_range.
         ax : matplotlib.axes._subplots.Axes = None
             Matplotlib axes to plot on, if not specified, one is created.
-        complex_fields : bool
-            Whether time domain fields are complex, e.g., for Bloch boundaries
 
         Returns
         -------
@@ -72,7 +70,7 @@ def plot_spectrum(
 
         fmin, fmax = self.frequency_range()
         return self.plot_spectrum_in_frequency_range(
-            times, fmin, fmax, num_freqs=num_freqs, val=val, ax=ax, complex_fields=complex_fields
+            times, fmin, fmax, num_freqs=num_freqs, val=val, ax=ax
         )
 
     @abstractmethod

tidy3d/components/source.py

@@ -40,17 +40,17 @@ def amp_time(self, time: float) -> complex:
         Returns
         -------
         complex
-            Complex-valued amplitude at that time..
+            Complex-valued amplitude at that time.
         """
 
     def spectrum(
         self,
         times: ArrayFloat1D,
         freqs: ArrayFloat1D,
         dt: float,
-        complex_fields: bool = False,
     ) -> complex:
-        """Complex-valued spectrum as a function of frequency
+        """Complex-valued spectrum as a function of frequency.
+        Note: Only the real part of the time signal is used.
 
         Parameters
         ----------
@@ -62,8 +62,6 @@ def spectrum(
         dt : float or np.ndarray
             Time step to weight FT integral with.
             If array, use to weigh each of the time intervals in ``times``.
-        complex_fields : bool
-            Whether time domain fields are complex, e.g., for Bloch boundaries
 
         Returns
         -------
@@ -73,10 +71,7 @@ def spectrum(
 
         times = np.array(times)
         freqs = np.array(freqs)
-        time_amps = self.amp_time(times)
-
-        if not complex_fields:
-            time_amps = np.real(time_amps)
+        time_amps = np.real(self.amp_time(times))
 
         # if all time amplitudes are zero, just return (complex-valued) zeros for spectrum
         if np.all(np.equal(time_amps, 0.0)):
@@ -86,7 +81,7 @@ def spectrum(
         relevant_time_inds = np.where(np.abs(time_amps) / np.amax(np.abs(time_amps)) > DFT_CUTOFF)
         # find first and last index where the filter is True
         start_ind = relevant_time_inds[0][0]
-        stop_ind = relevant_time_inds[0][-1]
+        stop_ind = relevant_time_inds[0][-1] + 1
         time_amps = time_amps[start_ind:stop_ind]
         times_cut = times[start_ind:stop_ind]
         if times_cut.size == 0:
@@ -118,9 +113,9 @@ def plot_spectrum_in_frequency_range(
         num_freqs: int = 101,
         val: PlotVal = "real",
         ax: Ax = None,
-        complex_fields: bool = False,
     ) -> Ax:
         """Plot the complex-valued amplitude of the time-dependence.
+        Note: Only the real part of the time signal is used.
 
         Parameters
         ----------
@@ -137,8 +132,6 @@ def plot_spectrum_in_frequency_range(
             Number of frequencies to plot within the [fmin, fmax].
         ax : matplotlib.axes._subplots.Axes = None
             Matplotlib axes to plot on, if not specified, one is created.
-        complex_fields : bool
-            Whether time domain fields are complex, e.g., for Bloch boundaries
 
         Returns
         -------
@@ -154,7 +147,7 @@ def plot_spectrum_in_frequency_range(
         dt = np.mean(dts)
         freqs = np.linspace(fmin, fmax, num_freqs)
 
-        spectrum = self.spectrum(times=times, dt=dt, freqs=freqs, complex_fields=complex_fields)
+        spectrum = self.spectrum(times=times, dt=dt, freqs=freqs)
 
         if val == "real":
             ax.plot(freqs, spectrum.real, color="blueviolet", label="real")