diff --git a/CHANGELOG.md b/CHANGELOG.md index 5b6a541e3..8d076274f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Renamed `SkinDepthFitterParam` --> `SurfaceImpedanceFitterParam` used in `LossyMetalMedium`. `plot` in `LossyMetalMedium` now plots complex-valued surface impedance. - Changed plot_3d iframe url to tidy3d production environment. - `num_freqs` is now set to 3 by default for the `PlaneWave`, `GaussianBeam`, and `AnalyticGaussianBeam` sources, which makes the injection more accurate in broadband cases. +- Nonlinear models `KerrNonlinearity` and `TwoPhotonAbsorption` now default to using the physical real fields instead of complex fields. ### Fixed - Make gauge selection for non-converged modes more robust. diff --git a/tests/test_components/test_medium.py b/tests/test_components/test_medium.py index b72548f63..ec6cae163 100644 --- a/tests/test_components/test_medium.py +++ b/tests/test_components/test_medium.py @@ -591,6 +591,18 @@ def test_nonlinear_medium(): ) ) + assert med._nonlinear_num_iters == 20 + assert td.Medium()._nonlinear_num_iters == 0 + assert td.Medium(nonlinear_spec=td.NonlinearSusceptibility(chi3=1.5))._nonlinear_num_iters == 1 + assert ( + td.Medium( + nonlinear_spec=td.NonlinearSusceptibility(chi3=1.5, numiters=2) + )._nonlinear_num_iters + == 2 + ) + assert td.Medium()._nonlinear_models == [] + assert td.Medium(nonlinear_spec=td.NonlinearSpec())._nonlinear_models == [] + # warn about deprecated api with AssertLogLevel("WARNING"): med = td.Medium(nonlinear_spec=td.NonlinearSusceptibility(chi3=1.5)) @@ -634,11 +646,15 @@ def test_nonlinear_medium(): # active materials with pytest.raises(ValidationError): med = td.Medium( - nonlinear_spec=td.NonlinearSpec(models=[td.TwoPhotonAbsorption(beta=-1, n0=1)]) + nonlinear_spec=td.NonlinearSpec(models=[td.TwoPhotonAbsorption(beta=-1, n0=1, freq0=1)]) ) with pytest.raises(ValidationError): - med = td.Medium(nonlinear_spec=td.NonlinearSpec(models=[td.KerrNonlinearity(n2=-1j, n0=1)])) + med = td.Medium( + nonlinear_spec=td.NonlinearSpec( + models=[td.KerrNonlinearity(n2=-1j, n0=1, use_complex_fields=True)] + ) + ) # automatic detection of n0 and freq0 n0 = 2 @@ -662,26 +678,6 @@ def test_nonlinear_medium(): assert n0 == nonlinear_spec.models[0]._get_n0(n0=None, medium=medium, freqs=[freq0]) assert freq0 == nonlinear_spec.models[0]._get_freq0(freq0=None, freqs=[freq0]) - # two photon absorption is phenomenological - with AssertLogLevel("WARNING", contains_str="phenomenological"): - med = td.Medium( - nonlinear_spec=td.NonlinearSpec(models=[td.TwoPhotonAbsorption(beta=-1, n0=1)]), - allow_gain=True, - ) - sim.updated_copy(medium=med, path="structures/0") - - # complex parameters - with AssertLogLevel("WARNING", contains_str="preferred"): - med = td.Medium( - nonlinear_spec=td.NonlinearSpec( - models=[ - td.KerrNonlinearity(n2=-1 + 1j, n0=1), - ], - num_iters=20, - ) - ) - sim.updated_copy(medium=med, path="structures/0") - # subsection with nonlinear materials needs to hardcode source info sim2 = sim.updated_copy(center=(-4, -4, -4), path="sources/0") sim2 = sim2.updated_copy( @@ -696,6 +692,8 @@ def test_nonlinear_medium(): source2 = source.updated_copy(source_time=source_time2) with pytest.raises(SetupError): sim.updated_copy(sources=[source, source2]) + with pytest.raises(SetupError): + sim.updated_copy(sources=[]) # but if we provided it, it's ok nonlinear_spec = td.NonlinearSpec(models=[td.KerrNonlinearity(n2=1, n0=1)]) @@ -703,6 +701,18 @@ def test_nonlinear_medium(): sim = sim.updated_copy(structures=[structure]) assert 1 == nonlinear_spec.models[0]._get_n0(n0=1, medium=medium, freqs=[1, 2]) + nonlinear_spec = td.NonlinearSpec(models=[td.TwoPhotonAbsorption(beta=1, n0=1)]) + structure = structure.updated_copy(medium=medium.updated_copy(nonlinear_spec=nonlinear_spec)) + sim = sim.updated_copy(structures=[structure]) + with pytest.raises(SetupError): + sim = sim.updated_copy(structures=[structure], sources=[source, source2]) + with pytest.raises(SetupError): + sim = sim.updated_copy(structures=[structure], sources=[]) + nonlinear_spec = td.NonlinearSpec(models=[td.TwoPhotonAbsorption(beta=1, n0=1, freq0=1)]) + structure = structure.updated_copy(medium=medium.updated_copy(nonlinear_spec=nonlinear_spec)) + sim = sim.updated_copy(structures=[structure]) + assert 1 == nonlinear_spec.models[0]._get_freq0(freq0=1, freqs=[1, 2]) + # active materials with automatic detection of n0 nonlinear_spec_active = td.NonlinearSpec(models=[td.TwoPhotonAbsorption(beta=-1)]) medium_active = medium.updated_copy(nonlinear_spec=nonlinear_spec_active) @@ -727,6 +737,38 @@ def test_nonlinear_medium(): with pytest.raises(ValidationError): td.Medium2D(ss=modulated, tt=modulated) + # some parameters must be real now, unless we use old implementation + _ = td.TwoPhotonAbsorption(beta=1j, use_complex_fields=True) + with pytest.raises(pydantic.ValidationError): + _ = td.TwoPhotonAbsorption(beta=1j) + _ = td.KerrNonlinearity(n2=1j, use_complex_fields=True) + with pytest.raises(pydantic.ValidationError): + _ = td.KerrNonlinearity(n2=1j) + + # consistent complex fields + _ = td.NonlinearSpec( + models=[ + td.TwoPhotonAbsorption(beta=1, use_complex_fields=True), + td.KerrNonlinearity(n2=1, use_complex_fields=True), + ] + ) + with pytest.raises(pydantic.ValidationError): + _ = td.NonlinearSpec( + models=[ + td.TwoPhotonAbsorption(beta=1, use_complex_fields=True), + td.KerrNonlinearity(n2=1, use_complex_fields=False), + ] + ) + + # warn if using old implementation + with AssertLogLevel("WARNING", contains_str="use_complex_fields"): + med = td.Medium( + nonlinear_spec=td.NonlinearSpec( + models=[td.KerrNonlinearity(n2=1, use_complex_fields=True)] + ) + ) + _ = sim.updated_copy(medium=med, path="structures/0") + def test_custom_medium(): Nx, Ny, Nz, Nf = 4, 3, 1, 1 diff --git a/tests/test_components/test_simulation.py b/tests/test_components/test_simulation.py index 97070ea3a..a4324720a 100644 --- a/tests/test_components/test_simulation.py +++ b/tests/test_components/test_simulation.py @@ -3073,7 +3073,7 @@ def test_fixed_angle_sim(): permittivity=3, nonlinear_spec=td.NonlinearSpec( models=[ - td.KerrNonlinearity(n2=-1 + 1j, n0=1), + td.KerrNonlinearity(n2=1, n0=1), ], num_iters=20, ), diff --git a/tidy3d/components/medium.py b/tidy3d/components/medium.py index cd5cae09d..0dacc3c4f 100644 --- a/tidy3d/components/medium.py +++ b/tidy3d/components/medium.py @@ -218,6 +218,8 @@ def _get_n0( freqs: List[pd.PositiveFloat], ) -> complex: """Get a single value for n0.""" + if freqs is None: + freqs = [] freqs = np.array(freqs, dtype=float) ns, ks = medium.nk_model(freqs) nks = ns + 1j * ks @@ -256,7 +258,7 @@ def _get_n0( @property def complex_fields(self) -> bool: """Whether the model uses complex fields.""" - pass + return False class NonlinearSusceptibility(NonlinearModel): @@ -322,11 +324,6 @@ 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 @@ -334,41 +331,60 @@ class TwoPhotonAbsorption(NonlinearModel): Also includes free-carrier absorption (FCA) and free-carrier plasma dispersion (FCPD) effects. The expression for the nonlinear polarization is given below. - Note - ---- - .. math:: + Notes + ----- - P_{NL} = P_{TPA} + P_{FCA} + P_{FCPD} \\\\ - P_{TPA} = -\\frac{c_0^2 \\varepsilon_0^2 n_0 \\operatorname{Re}(n_0) \\beta}{2 i \\omega} |E|^2 E \\\\ - P_{FCA} = -\\frac{c_0 \\varepsilon_0 n_0 \\sigma N_f}{i \\omega} E \\\\ - \\frac{dN_f}{dt} = \\frac{c_0^2 \\varepsilon_0^2 n_0^2 \\beta}{8 q_e \\hbar \\omega} |E|^4 - \\frac{N_f}{\\tau} \\\\ - N_e = N_h = N_f \\\\ - P_{FCPD} = \\varepsilon_0 2 n_0 \\Delta n (N_f) E \\\\ - \\Delta n (N_f) = (c_e N_e^{e_e} + c_h N_h^{e_h}) + This model uses real time-domain fields, so :math:`\\beta` must be real. - Note - ---- - This frequency-domain equation is implemented in the time domain using complex-valued fields. + .. math:: - Note - ---- - Different field components do not interact nonlinearly. For example, - when calculating :math:`P_{NL, x}`, we approximate :math:`|E|^2 \\approx |E_x|^2`. - This approximation is valid when the E field is predominantly polarized along one - of the x, y, or z axes. + P_{NL} = P_{TPA} + P_{FCA} + P_{FCPD} \\\\ + P_{TPA} = -\\frac{4}{3}\\frac{c_0^2 \\varepsilon_0^2 n_0^2 \\beta}{2 i \\omega} |E|^2 E \\\\ + P_{FCA} = -\\frac{c_0 \\varepsilon_0 n_0 \\sigma N_f}{i \\omega} E \\\\ + \\frac{dN_f}{dt} = \\frac{8}{3}\\frac{c_0^2 \\varepsilon_0^2 n_0^2 \\beta}{8 q_e \\hbar \\omega} |E|^4 - \\frac{N_f}{\\tau} \\\\ + N_e = N_h = N_f \\\\ + P_{FCPD} = \\varepsilon_0 2 n_0 \\Delta n (N_f) E \\\\ + \\Delta n (N_f) = (c_e N_e^{e_e} + c_h N_h^{e_h}) - Note - ---- - The implementation is described in:: + In these equations, :math:`n_0` means the real part of the linear + refractive index of the medium. + + The nonlinear constitutive relation is solved iteratively; it may not converge + for strong nonlinearities. Increasing :attr:`tidy3d.NonlinearSpec.num_iters` can + help with convergence. + + For complex fields (e.g. when using Bloch boundary conditions), the nonlinearity + is applied separately to the real and imaginary parts, so that the above equation + holds when both :math:`E` and :math:`P_{NL}` are replaced by their real or imaginary parts. + The nonlinearity is only applied to the real-valued fields since they are the + physical fields. + + Different field components do not interact nonlinearly. For example, + when calculating :math:`P_{NL, x}`, we approximate :math:`|E|^2 \\approx |E_x|^2`. + This approximation is valid when the :math:`E` field is predominantly polarized along one + of the ``x``, ``y``, or ``z`` axes. - N. Suzuki, "FDTD Analysis of Two-Photon Absorption and Free-Carrier Absorption in Si - High-Index-Contrast Waveguides," J. Light. Technol. 25, 9 (2007). + The implementation is described in:: + + N. Suzuki, "FDTD Analysis of Two-Photon Absorption and Free-Carrier Absorption in Si + High-Index-Contrast Waveguides," J. Light. Technol. 25, 9 (2007). + + .. TODO add links to notebooks here. Example ------- >>> tpa_model = TwoPhotonAbsorption(beta=1) """ + use_complex_fields: bool = pd.Field( + False, + title="Use complex fields", + description="Whether to use the old deprecated complex-fields implementation. " + "The default real-field implementation is more physical and is always " + "recommended; this option is only available for backwards compatibility " + "with Tidy3D version < 2.8 and may be removed in a future release.", + ) + beta: Union[float, Complex] = pd.Field( 0, title="TPA coefficient", @@ -429,10 +445,26 @@ class TwoPhotonAbsorption(NonlinearModel): "from the simulation sources (as long as these are all equal).", ) + @pd.validator("beta", always=True) + def _validate_beta_real(cls, val, values): + """Check that beta is real and give a useful error if it is not.""" + use_complex_fields = values.get("use_complex_fields") + if use_complex_fields: + return val + if not np.isreal(val): + raise SetupError( + "Complex values of 'beta' in 'TwoPhotonAbsorption' are not " + "supported; the implementation uses the " + "physical real-valued fields." + ) + return val + def _validate_medium_freqs(self, medium: AbstractMedium, freqs: List[pd.PositiveFloat]) -> None: """Any validation that depends on knowing the central frequencies of the sources. This includes passivity checking, if necessary.""" n0 = self._get_n0(self.n0, medium, freqs) + if freqs is not None: + _ = self._get_freq0(self.freq0, freqs) beta = self.beta if not medium.allow_gain: chi_imag = np.real(beta * n0 * np.real(n0)) @@ -457,12 +489,12 @@ def _validate_medium(self, medium: AbstractMedium): """Check that the model is compatible with the medium.""" # if n0 is specified, we can go ahead and validate passivity if self.n0 is not None: - self._validate_medium_freqs(medium, []) + self._validate_medium_freqs(medium, None) @property def complex_fields(self) -> bool: """Whether the model uses complex fields.""" - return True + return self.use_complex_fields class KerrNonlinearity(NonlinearModel): @@ -470,34 +502,57 @@ class KerrNonlinearity(NonlinearModel): of the form :math:`n = n_0 + n_2 I`. The expression for the nonlinear polarization is given below. - Note - ---- - .. math:: + Notes + ----- - P_{NL} = \\varepsilon_0 c_0 n_0 \\operatorname{Re}(n_0) n_2 |E|^2 E + This model uses real time-domain fields, so :math:`\\n_2` must be real. - Note - ---- - The fields in this equation are complex-valued, allowing a direct implementation of the Kerr - nonlinearity. In contrast, the model :class:`.NonlinearSusceptibility` implements a - chi3 nonlinear susceptibility using real-valued fields, giving rise to Kerr nonlinearity - as well as third-harmonic generation and other effects. The relationship between the parameters is given by - :math:`n_2 = \\frac{3}{4} \\frac{1}{\\varepsilon_0 c_0 n_0 \\operatorname{Re}(n_0)} \\chi_3`. The additional - factor of :math:`\\frac{3}{4}` comes from the usage of complex-valued fields for the Kerr - nonlinearity and real-valued fields for the nonlinear susceptibility. + This model is equivalent to a :class:`.NonlinearSusceptibility`; the + relation between the parameters is given below. - Note - ---- - Different field components do not interact nonlinearly. For example, - when calculating :math:`P_{NL, x}`, we approximate :math:`|E|^2 \\approx |E_x|^2`. - This approximation is valid when the E field is predominantly polarized along one - of the x, y, or z axes. + .. math:: + + P_{NL} = \\varepsilon_0 \\chi_3 |E|^2 E \\\\ + n_2 = \\frac{3}{4 n_0^2 \\varepsilon_0 c_0} \\chi_3 + + In these equations, :math:`n_0` means the real part of the linear + refractive index of the medium. + + To simulate nonlinear loss, consider instead using a :class:`.TwoPhotonAbsorption` + model, which implements a more physical dispersive loss of the form + :math:`\\chi_{TPA} = i \\frac{c_0 n_0 \\beta}{\\omega} I`. + + The nonlinear constitutive relation is solved iteratively; it may not converge + for strong nonlinearities. Increasing :attr:`tidy3d.NonlinearSpec.num_iters` can + help with convergence. + + For complex fields (e.g. when using Bloch boundary conditions), the nonlinearity + is applied separately to the real and imaginary parts, so that the above equation + holds when both :math:`E` and :math:`P_{NL}` are replaced by their real or imaginary parts. + The nonlinearity is only applied to the real-valued fields since they are the + physical fields. + + Different field components do not interact nonlinearly. For example, + when calculating :math:`P_{NL, x}`, we approximate :math:`|E|^2 \\approx |E_x|^2`. + This approximation is valid when the :math:`E` field is predominantly polarized along one + of the ``x``, ``y``, or ``z`` axes. + + .. TODO add links to notebooks here. Example ------- >>> kerr_model = KerrNonlinearity(n2=1) """ + use_complex_fields: bool = pd.Field( + False, + title="Use complex fields", + description="Whether to use the old deprecated complex-fields implementation. " + "The default real-field implementation is more physical and is always " + "recommended; this option is only available for backwards compatibility " + "with Tidy3D version < 2.8 and may be removed in a future release.", + ) + n2: Complex = pd.Field( 0, title="Nonlinear refractive index", @@ -513,11 +568,31 @@ class KerrNonlinearity(NonlinearModel): "frequencies of the simulation sources (as long as these are all equal).", ) + @pd.validator("n2", always=True) + def _validate_n2_real(cls, val, values): + """Check that n2 is real and give a useful error if it is not.""" + use_complex_fields = values.get("use_complex_fields") + if use_complex_fields: + return val + if not np.isreal(val): + raise SetupError( + "Complex values of 'n2' in 'KerrNonlinearity' are not " + "supported; the implementation uses the " + "physical real-valued fields. " + "To simulate nonlinear loss, consider instead using a " + "'TwoPhotonAbsorption' model, which implements a " + "more physical dispersive loss of the form " + "'chi_{TPA} = i (c_0 n_0 beta / omega) I'." + ) + return val + def _validate_medium_freqs(self, medium: AbstractMedium, freqs: List[pd.PositiveFloat]) -> None: """Any validation that depends on knowing the central frequencies of the sources. This includes passivity checking, if necessary.""" n0 = self._get_n0(self.n0, medium, freqs) n2 = self.n2 + if not self.use_complex_fields: + return if not medium.allow_gain: chi_imag = np.imag(n2 * n0 * np.real(n0)) if chi_imag < 0: @@ -529,6 +604,12 @@ def _validate_medium_freqs(self, medium: AbstractMedium, freqs: List[pd.Positive "gain medium are unstable, and are likely to diverge." ) + def _validate_medium(self, medium: AbstractMedium): + """Check that the model is compatible with the medium.""" + # if n0 is specified, we can go ahead and validate passivity + if self.n0 is not None: + self._validate_medium_freqs(medium, []) + def _hardcode_medium_freqs( self, medium: AbstractMedium, freqs: List[pd.PositiveFloat] ) -> KerrNonlinearity: @@ -536,16 +617,10 @@ def _hardcode_medium_freqs( n0 = self._get_n0(n0=self.n0, medium=medium, freqs=freqs) return self.updated_copy(n0=n0) - def _validate_medium(self, medium: AbstractMedium): - """Check that the model is compatible with the medium.""" - # if n0 is specified, we can go ahead and validate passivity - 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 + return self.use_complex_fields NonlinearModelType = Union[NonlinearSusceptibility, TwoPhotonAbsorption, KerrNonlinearity] @@ -595,6 +670,29 @@ def _no_duplicate_models(cls, val): ) return val + @pd.validator("models", always=True) + def _consistent_old_complex_fields(cls, val): + """Ensure that old complex fields implementation is used consistently.""" + if val is None: + return val + use_complex_fields = False + for model in val: + if isinstance(model, (KerrNonlinearity, TwoPhotonAbsorption)): + if model.use_complex_fields: + use_complex_fields = True + elif use_complex_fields: + # if one model uses complex fields, they all should + raise SetupError( + "Some of the nonlinear models have " + "'use_complex_fields=True' and some have " + "'use_complex_fields=False'. This option " + "is only available for backwards compatibility " + "with Tidy3D version < 2.8 " + "and it must be consistent across the nonlinear " + "models in a given 'NonlinearSpec'." + ) + return val + @pd.validator("num_iters", always=True) def _validate_num_iters(cls, val, values): """Check that num_iters is not too large.""" @@ -657,7 +755,7 @@ class AbstractMedium(ABC, Tidy3dBaseModel): ) @cached_property - def _nonlinear_models(self) -> NonlinearSpec: + def _nonlinear_models(self) -> List: """The nonlinear models in the nonlinear_spec.""" if self.nonlinear_spec is None: return [] @@ -665,7 +763,7 @@ def _nonlinear_models(self) -> NonlinearSpec: return [self.nonlinear_spec] if self.nonlinear_spec.models is None: return [] - return self.nonlinear_spec.models + return list(self.nonlinear_spec.models) @cached_property def _nonlinear_num_iters(self) -> pd.PositiveInt: diff --git a/tidy3d/components/simulation.py b/tidy3d/components/simulation.py index 9a64abed7..c55f2e849 100644 --- a/tidy3d/components/simulation.py +++ b/tidy3d/components/simulation.py @@ -57,13 +57,11 @@ AbstractPerturbationMedium, AnisotropicMedium, FullyAnisotropicMedium, - KerrNonlinearity, LossyMetalMedium, Medium, Medium2D, MediumType, MediumType3D, - TwoPhotonAbsorption, ) from .monitor import ( AbstractFieldProjectionMonitor, @@ -3593,30 +3591,12 @@ def _validate_nonlinear_specs(self) -> None: for model in medium._nonlinear_models: model._validate_medium_freqs(medium, freqs) - if isinstance(model, TwoPhotonAbsorption): - if np.iscomplex(model.beta): - log.warning( - "Complex values of 'beta' in 'TwoPhotonAbsorption' are deprecated " - "and may be removed in a future version. The implementation with " - "complex 'beta' is as described in the 'TwoPhotonAbsorption' docstring, " - "but the physical interpretation of 'beta' may not be correct if it is complex." - ) - log.warning( - "Found a medium with a 'TwoPhotonAbsorption' nonlinearity. " - "This uses a phenomenological model based on complex fields, " - "so care should be taken in interpreting the results. For more " - "information on the model, see the documentation at " - "'https://docs.flexcompute.com/projects/tidy3d/en/latest/api/_autosummary/tidy3d.TwoPhotonAbsorption.html' or the following reference: " - "'N. Suzuki, \"FDTD Analysis of Two-Photon Absorption and Free-Carrier Absorption in Si High-Index-Contrast Waveguides,\" J. Light. Technol. 25, 9 (2007).'." - ) - - if isinstance(model, KerrNonlinearity): + if model.complex_fields: log.warning( - "Found a medium with a 'KerrNonlinearity'. Usually, " - "'NonlinearSusceptibility' is preferred, as it captures " - "additional physical effects by acting on the underlying real fields. " - "The relation between the parameters is " - "'chi3 = (4/3) * eps_0 * c_0 * n0 * Re(n0) * n2'." + "Found a nonlinear model with 'use_complex_fields=True'. " + "For physical simulation results, this should always " + "be 'False'. This option is available only for backwards " + "compatibility and may be removed in a future release." ) """ Pre submit validation (before web.upload()) """