From 63668c78db5b6d7a06b15966b4cbb8c147bb39c4 Mon Sep 17 00:00:00 2001 From: Casey Wojcik Date: Wed, 6 Sep 2023 04:35:37 +0100 Subject: [PATCH] Added TwoPhotonAbsorption and KerrNonlinearity --- CHANGELOG.md | 3 + tests/sims/simulation_2_5_0rc2.json | 58 +++- tests/test_components/test_medium.py | 119 +++++++- tests/test_components/test_source.py | 35 ++- tests/utils.py | 16 ++ tidy3d/__init__.py | 7 +- tidy3d/components/medium.py | 399 +++++++++++++++++++++++---- tidy3d/components/simulation.py | 27 +- tidy3d/components/source.py | 12 + tidy3d/constants.py | 2 + 10 files changed, 597 insertions(+), 81 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 4e60db699..b3e851eb2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,9 +6,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] ### Added +- Added support for two-photon absorption via `TwoPhotonAbsorption` class. Added `KerrNonlinearity` that implements Kerr effect without third-harmonic generation. ### Changed - Indent for the json string of Tidy3D models has been changed to `None` when used internally; kept as `indent=4` for writing to `json` and `yaml` files. +- API for specifying one or more nonlinear models via `NonlinearSpec.models`. ### Fixed - Fixed the duplication of log messages in Jupyter when `set_logging_file` is used. @@ -32,6 +34,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - python 3.7 no longer tested nor supported. - Removed warning that monitors now have `colocate=True` by default. - If `PML` or any absorbing boundary condition is used along a direction where the `Simulation` size is zero, an error will be raised, rather than just a warning. +- Remove warning that monitors now have `colocate=True` by default. ### Fixed - If there are no adjoint sources for a simulation involved in an objective function, make a mock source with zero amplitude and warn user. diff --git a/tests/sims/simulation_2_5_0rc2.json b/tests/sims/simulation_2_5_0rc2.json index 70702aab5..127b7761a 100644 --- a/tests/sims/simulation_2_5_0rc2.json +++ b/tests/sims/simulation_2_5_0rc2.json @@ -675,9 +675,63 @@ "frequency_range": null, "allow_gain": false, "nonlinear_spec": { - "numiters": 20, "type": "NonlinearSusceptibility", - "chi3": 0.1 + "chi3": 0.1, + "numiters": 20 + }, + "modulation_spec": null, + "heat_spec": null, + "type": "Medium", + "permittivity": 1.0, + "conductivity": 0.0 + } + }, + { + "geometry": { + "type": "Box", + "center": [ + -1.0, + 0.5, + 0.5 + ], + "size": [ + 1.0, + 1.0, + 1.0 + ] + }, + "name": null, + "type": "Structure", + "medium": { + "name": null, + "frequency_range": null, + "allow_gain": false, + "nonlinear_spec": { + "models": [ + { + "type": "NonlinearSusceptibility", + "chi3": 0.1, + "numiters": null + }, + { + "type": "TwoPhotonAbsorption", + "beta": { + "real": 1.0, + "imag": 0.0 + }, + "n0": null + }, + { + "type": "KerrNonlinearity", + "n2": { + "real": 1.0, + "imag": 0.0 + }, + "n0": null + } + ], + "num_iters": 10, + "type": "NonlinearSpec" }, "modulation_spec": null, "heat_spec": null, diff --git a/tests/test_components/test_medium.py b/tests/test_components/test_medium.py index e8b74f0a7..b0d561bb0 100644 --- a/tests/test_components/test_medium.py +++ b/tests/test_components/test_medium.py @@ -4,7 +4,7 @@ import pydantic.v1 as pydantic import matplotlib.pyplot as plt import tidy3d as td -from tidy3d.exceptions import ValidationError +from tidy3d.exceptions import ValidationError, SetupError from ..utils import assert_log_level, log_capture from typing import Dict @@ -546,13 +546,120 @@ def test_perturbation_medium(): ) -def test_nonlinear_medium(): - med = td.Medium(nonlinear_spec=td.NonlinearSusceptibility(chi3=1.5, numiters=20)) +def test_nonlinear_medium(log_capture): + med = td.Medium( + nonlinear_spec=td.NonlinearSpec( + models=[ + td.NonlinearSusceptibility(chi3=1.5), + td.TwoPhotonAbsorption(beta=1), + td.KerrNonlinearity(n2=1), + ], + num_iters=20, + ) + ) - with pytest.raises(pydantic.ValidationError): - med = td.PoleResidue( - poles=[(-1, 1)], nonlinear_spec=td.NonlinearSusceptibility(chi3=1.5, numiters=20) + # complex parameters + med = td.Medium( + nonlinear_spec=td.NonlinearSpec( + models=[ + td.KerrNonlinearity(n2=-1 + 1j, n0=1), + td.TwoPhotonAbsorption(beta=1 + 1j, n0=1), + ], + num_iters=20, + ) + ) + assert_log_level(log_capture, None) + + # warn about deprecated api + med = td.Medium(nonlinear_spec=td.NonlinearSusceptibility(chi3=1.5)) + assert_log_level(log_capture, "WARNING") + + # don't use deprecated numiters + with pytest.raises(ValidationError): + med = td.Medium( + nonlinear_spec=td.NonlinearSpec(models=[td.NonlinearSusceptibility(chi3=1, numiters=2)]) + ) + + # dispersive support + med = td.PoleResidue(poles=[(-1, 1)], nonlinear_spec=td.NonlinearSusceptibility(chi3=1.5)) + + # unsupported material types + with pytest.raises(ValidationError): + med = td.AnisotropicMedium( + xx=med, yy=med, zz=med, nonlinear_spec=td.NonlinearSusceptibility(chi3=1.5) ) + # numiters too large with pytest.raises(pydantic.ValidationError): med = td.Medium(nonlinear_spec=td.NonlinearSusceptibility(chi3=1.5, numiters=200)) + with pytest.raises(pydantic.ValidationError): + med = td.Medium( + nonlinear_spec=td.NonlinearSpec( + num_iters=200, models=[td.NonlinearSusceptibility(chi3=1.5)] + ) + ) + + # duplicate models + with pytest.raises(pydantic.ValidationError): + med = td.Medium( + nonlinear_spec=td.NonlinearSpec( + models=[ + td.NonlinearSusceptibility(chi3=1.5), + td.NonlinearSusceptibility(chi3=1), + ] + ) + ) + + # active materials + with pytest.raises(ValidationError): + med = td.Medium( + nonlinear_spec=td.NonlinearSpec(models=[td.TwoPhotonAbsorption(beta=-1 + 1j, n0=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.TwoPhotonAbsorption(beta=-1, n0=1)]), + allow_gain=True, + ) + + # automatic detection of n0 + n0 = 2 + freq0 = td.C_0 / 1 + nonlinear_spec = td.NonlinearSpec(models=[td.KerrNonlinearity(n2=1)]) + medium = td.Sellmeier.from_dispersion(n=n0, freq=freq0, dn_dwvl=-0.2).updated_copy( + nonlinear_spec=nonlinear_spec + ) + source_time = td.GaussianPulse(freq0=freq0, fwidth=freq0 / 10) + source = td.PointDipole(center=(0, 0, 0), source_time=source_time, polarization="Ex") + monitor = td.FieldMonitor(size=(td.inf, td.inf, 0), freqs=[freq0], name="field") + structure = td.Structure(geometry=td.Box(size=(5, 5, 5)), medium=medium) + sim = td.Simulation( + size=(10, 10, 10), + run_time=1e-12, + grid_spec=td.GridSpec.uniform(dl=0.1), + sources=[source], + monitors=[monitor], + structures=[structure], + ) + assert n0 == nonlinear_spec.models[0]._get_n0(n0=None, medium=medium, freqs=[freq0]) + + # can't detect n0 with different source freqs + source_time2 = source_time.updated_copy(freq0=2 * freq0) + source2 = source.updated_copy(source_time=source_time2) + with pytest.raises(SetupError): + sim.updated_copy(sources=[source, source2]) + + # but if we provided it, it's ok + nonlinear_spec = td.NonlinearSpec(models=[td.KerrNonlinearity(n2=1, n0=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_n0(n0=1, medium=medium, 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) + with pytest.raises(ValidationError): + structure = structure.updated_copy(medium=medium_active) + sim.updated_copy(structures=[structure]) diff --git a/tests/test_components/test_source.py b/tests/test_components/test_source.py index c59a5868a..a756df266 100644 --- a/tests/test_components/test_source.py +++ b/tests/test_components/test_source.py @@ -6,6 +6,7 @@ import tidy3d as td from tidy3d.exceptions import SetupError from tidy3d.components.source import DirectionalSource, CHEB_GRID_WIDTH +from ..utils import assert_log_level, log_capture ST = td.GaussianPulse(freq0=2e14, fwidth=1e14) S = td.PointDipole(source_time=ST, polarization="Ex") @@ -268,7 +269,7 @@ def check_freq_grid(freq_grid, num_freqs): ) -def test_custom_source_time(): +def test_custom_source_time(log_capture): ts = np.linspace(0, 30, 1001) amp_time = ts / max(ts) @@ -276,31 +277,22 @@ def test_custom_source_time(): 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 * np.exp(-1j * 2 * np.pi * ts), 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 * np.exp(-1j * 2 * np.pi * 0.09)], 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) + # test out of range handling source = td.PointDipole(center=(0, 0, 0), source_time=cst, polarization="Ex") + monitor = td.FieldMonitor(size=(td.inf, td.inf, 0), freqs=[1], name="field") sim = td.Simulation( size=(10, 10, 10), run_time=1e-12, grid_spec=td.GridSpec.uniform(dl=0.1), sources=[source], + normalize_index=None, ) - - # test out of range handling - vals = [1] cst = td.CustomSourceTime.from_values(freq0=1, fwidth=0.1, values=[0, 1], dt=sim.dt) source = td.PointDipole(center=(0, 0, 0), source_time=cst, polarization="Ex") sim = sim.updated_copy(sources=[source]) @@ -308,3 +300,20 @@ def test_custom_source_time(): assert np.allclose( cst.amp_time(sim.tmesh[1:]), np.exp(-1j * 2 * np.pi * sim.tmesh[1:]), rtol=0, atol=ATOL ) + + assert_log_level(log_capture, None) + + # test normalization warning + sim = sim.updated_copy(normalize_index=0) + assert_log_level(log_capture, "WARNING") + log_capture.clear() + source = source.updated_copy(source_time=td.ContinuousWave(freq0=1, fwidth=0.1)) + sim = sim.updated_copy(sources=[source]) + assert_log_level(log_capture, "WARNING") + + # 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) diff --git a/tests/utils.py b/tests/utils.py index 01bf66aa5..3f10be28f 100644 --- a/tests/utils.py +++ b/tests/utils.py @@ -213,6 +213,22 @@ nonlinear_spec=td.NonlinearSusceptibility(chi3=0.1, numiters=20), ), ), + td.Structure( + geometry=td.Box( + size=(1, 1, 1), + center=(-1.0, 0.5, 0.5), + ), + medium=td.Medium( + nonlinear_spec=td.NonlinearSpec( + num_iters=10, + models=[ + td.NonlinearSusceptibility(chi3=0.1), + td.TwoPhotonAbsorption(beta=1), + td.KerrNonlinearity(n2=1), + ], + ) + ), + ), td.Structure( geometry=td.PolySlab( vertices=[(-1.5, -1.5), (-0.5, -1.5), (-0.5, -0.5)], slab_bounds=[-1, 1] diff --git a/tidy3d/__init__.py b/tidy3d/__init__.py index ce6f2f469..2aa463e01 100644 --- a/tidy3d/__init__.py +++ b/tidy3d/__init__.py @@ -16,7 +16,7 @@ from .components.medium import CustomMedium, CustomPoleResidue from .components.medium import CustomSellmeier, FullyAnisotropicMedium from .components.medium import CustomLorentz, CustomDrude, CustomDebye, CustomAnisotropicMedium -from .components.medium import NonlinearSusceptibility +from .components.medium import NonlinearSusceptibility, TwoPhotonAbsorption, KerrNonlinearity from .components.transformation import RotationAroundAxis from .components.medium import PerturbationMedium, PerturbationPoleResidue from .components.parameter_perturbation import ParameterPerturbation @@ -94,7 +94,7 @@ from .material_library.parametric_materials import Graphene # for docs -from .components.medium import AbstractMedium, NonlinearSpec +from .components.medium import AbstractMedium, NonlinearSpec, NonlinearModel from .components.geometry.base import Geometry from .components.source import Source, SourceTime from .components.monitor import Monitor @@ -187,7 +187,10 @@ def set_logging_level(level: str) -> None: "LinearChargePerturbation", "CustomChargePerturbation", "NonlinearSpec", + "NonlinearModel", "NonlinearSusceptibility", + "TwoPhotonAbsorption", + "KerrNonlinearity", "Structure", "MeshOverrideStructure", "ModeSpec", diff --git a/tidy3d/components/medium.py b/tidy3d/components/medium.py index f160bb1db..250dd3c28 100644 --- a/tidy3d/components/medium.py +++ b/tidy3d/components/medium.py @@ -14,7 +14,7 @@ from .grid.grid import Coords, Grid from .types import PoleAndResidue, Ax, FreqBound, TYPE_TAG_STR from .types import InterpMethod, Bound, ArrayComplex3D, ArrayFloat1D -from .types import Axis, TensorReal +from .types import Axis, TensorReal, Complex from .data.dataset import PermittivityDataset from .data.data_array import SpatialDataArray, ScalarFieldDataArray, DATA_ARRAY_MAP from .viz import add_ax_if_none @@ -22,6 +22,7 @@ from .validators import validate_name_str, validate_parameter_perturbation from ..constants import C_0, pec_val, EPSILON_0, LARGE_NUMBER, fp_eps, HBAR from ..constants import HERTZ, CONDUCTIVITY, PERMITTIVITY, RADPERSEC, MICROMETER, SECOND +from ..constants import WATT, VOLT from ..exceptions import ValidationError, SetupError from ..log import log from .transformation import RotationType @@ -36,7 +37,8 @@ FILL_VALUE = "extrapolate" # cap on number of nonlinear iterations -NONLINEAR_MAX_NUMITERS = 100 +NONLINEAR_MAX_NUM_ITERS = 100 +NONLINEAR_DEFAULT_NUM_ITERS = 5 # Range for checking upper bound of Im[eps], in addition to extrema method. # The range is in unit of eV and it's in log scale. @@ -51,7 +53,6 @@ def ensure_freq_in_range(eps_model: Callable[[float], complex]) -> Callable[[flo @functools.wraps(eps_model) def _eps_model(self, frequency: float) -> complex: """New eps_model function.""" - # evaluate infs and None as FREQ_EVAL_INF is_inf_scalar = isinstance(frequency, float) and np.isinf(frequency) if frequency is None or is_inf_scalar: @@ -83,54 +84,230 @@ def _eps_model(self, frequency: float) -> complex: """ Medium Definitions """ -class NonlinearSpec(ABC, Tidy3dBaseModel): - """Abstract specification for adding a nonlinearity to a medium. +class NonlinearModel(ABC, Tidy3dBaseModel): + """Abstract model for a nonlinear material response. + Used as part of a :class:`.NonlinearSpec`.""" + + def _validate_medium_type(self, medium: AbstractMedium): + """Check that the model is compatible with the medium.""" + if isinstance(medium, AbstractCustomMedium): + raise ValidationError( + f"'NonlinearModel' of class '{type(self).__name__}' is not currently supported " + f"for medium class '{type(medium).__name__}'." + ) + if medium.time_modulated: + raise ValidationError( + f"'NonlinearModel' of class '{type(self).__name__}' is not currently supported " + f"for time-modulated medium class '{type(medium).__name__}'." + ) + if not isinstance(medium, (Medium, DispersiveMedium)): + raise ValidationError( + f"'NonlinearModel' of class '{type(self).__name__}' is not currently supported " + f"for medium class '{type(medium).__name__}'." + ) + + def _validate_medium(self, medium: AbstractMedium): + """Any additional validation that depends on the medium""" + pass + + def _validate_medium_freqs(self, medium: AbstractMedium, freqs: List[pd.PositiveFloat]) -> None: + """Any additional validation that depends on the central frequencies of the sources.""" + pass + + def _get_n0( + self, + n0: complex, + medium: AbstractMedium, + freqs: List[pd.PositiveFloat], + ) -> complex: + """Get a single value for n0.""" + freqs = np.array(freqs, dtype=float) + ns, ks = medium.nk_model(freqs) + nks = ns + 1j * ks + + # n0 not specified; need to calculate it + if n0 is None: + if not len(nks): + raise SetupError( + f"Class '{type(self).__name__}' cannot determine 'n0' in the absence of " + "sources. Please either specify 'n0' or add sources to the simulation." + ) + if not all(np.isclose(nk, nks[0]) for nk in nks): + raise SetupError( + f"Class '{type(self).__name__}' cannot determine 'n0' because at the source " + f"frequencies '{freqs}' the complex refractive indices '{nks}' of the medium " + f"are not all equal. Please specify 'n0' in '{type(self).__name__}' " + "to match the complex refractive index of the medium at the desired " + "source central frequency." + ) + return nks[0] + + # now, n0 is specified; we use it, but warn if it might be inconsistent + if not all(np.isclose(nk, n0) for nk in nks): + log.warning( + f"Class '{type(self).__name__}' given 'n0={n0}'. At the source frequencies " + f"'{freqs}' the medium has complex refractive indices '{nks}'. In order " + "to obtain correct nonlinearity parameters, the provided refractive index " + "should agree with the complex refractive index at the source frequencies. " + "The provided value of 'n0' is being used; the resulting nonlinearity parameters " + "may be incorrect for those sources where the complex refractive index of the " + "medium is different from this value." + ) + return n0 + + +class NonlinearSusceptibility(NonlinearModel): + """Model for an instantaneous nonlinear chi3 susceptibility. + The expression for the instantaneous nonlinear polarization is given below. - Note + Note ---- - The nonlinear constitutive relation is solved iteratively; it may not converge - for strong nonlinearities. Increasing `numiters` can help with convergence. + .. math:: + + P_{NL} = \\varepsilon_0 \\chi_3 |E|^2 E + + Note + ---- + This model uses real time-domain fields, so :math:`\\chi_3` must be real. + 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 E and :math:`P_{NL}` are replaced by their real or imaginary parts. + The nonlinearity is applied to the real and imaginary components separately since + each of those represents a physical field. + + 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. + + Example + ------- + >>> nonlinear_susceptibility = NonlinearSusceptibility(chi3=1) """ + chi3: float = pd.Field( + 0, + title="Chi3", + description="Chi3 nonlinear susceptibility.", + units=f"{MICROMETER}^2 / {VOLT}^2", + ) + numiters: pd.PositiveInt = pd.Field( - 1, + None, title="Number of iterations", - description="Number of iterations for solving nonlinear constitutive relation.", + description="Deprecated. The old usage 'nonlinear_spec=model' with 'model.numiters' " + "is deprecated and will be removed in a future release. The new usage is " + r"'nonlinear_spec=NonlinearSpec(models=\[model], num_iters=num_iters)'. Under the new " + "usage, this parameter is ignored, and 'NonlinearSpec.num_iters' is used instead.", ) @pd.validator("numiters", always=True) def _validate_numiters(cls, val): """Check that numiters is not too large.""" - if val > NONLINEAR_MAX_NUMITERS: + if val is None: + return val + if val > NONLINEAR_MAX_NUM_ITERS: raise ValidationError( - "'NonlinearSpec.numiters' must be less than " - f"{NONLINEAR_MAX_NUMITERS}, currently {val}." + "'NonlinearSusceptibility.numiters' must be less than " + f"{NONLINEAR_MAX_NUM_ITERS}, currently {val}." ) return val -class NonlinearSusceptibility(NonlinearSpec): - """Specification adding an instantaneous nonlinear susceptibility to a medium. - The expression for the instantaneous nonlinear polarization is given below. +class TwoPhotonAbsorption(NonlinearModel): + """Model for two-photon absorption (TPA) nonlinearity which gives an intensity-dependent + absorption of the form :math:`\\alpha = \\alpha_0 + \\beta I`. + The expression for the nonlinear polarization is given below. Note ---- .. math:: - P_{NL} = \\epsilon_0 \\chi_3 |E|^2 E + P_{NL} = -\\frac{c_0^2 \\varepsilon_0^2 n_0 \\operatorname{Re}(n_0) \\beta}{2 i \\omega} |E|^2 E Note ---- - The nonlinear constitutive relation is solved iteratively; it may not converge - for strong nonlinearities. Increasing `numiters` can help with convergence. + This frequency-domain equation is implemented in the time domain using complex-valued fields. Note ---- - 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 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 E field is predominantly polarized along one + of the x, y, or z axes. + + Note + ---- + 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). + + Example + ------- + >>> tpa_model = TwoPhotonAbsorption(beta=1) + """ + + beta: Complex = pd.Field( + 0, + title="TPA coefficient", + description="Coefficient for two-photon absorption (TPA).", + units=f"{MICROMETER} / {WATT}", + ) + + n0: Optional[Complex] = pd.Field( + None, + title="Complex linear refractive index", + description="Complex linear refractive index of the medium, computed for instance using " + "'medium.nk_model'. If not provided, it is calculated automatically using the central " + "frequencies of the simulation sources (as long as these are all equal).", + ) + + 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) + beta = self.beta + if not medium.allow_gain: + chi_imag = np.real(beta * n0 * np.real(n0)) + if chi_imag < 0: + raise ValidationError( + "For passive medium, 'beta' in 'TwoPhotonAbsorption' must satisfy " + f"'Re(beta * n0 * Re(n0)) >= 0'. Currently, this quantity equals '{chi_imag}', " + f"and the linear index is 'n0={n0}'. To simulate gain medium, please set " + "'allow_gain=True' in the medium class. Caution: simulations containing " + "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, []) + + +class KerrNonlinearity(NonlinearModel): + """Model for Kerr nonlinearity which gives an intensity-dependent refractive index + of the form :math:`n = n_0 + n_2 I`. The expression for the nonlinear polarization + is given below. + + Note + ---- + .. math:: + + P_{NL} = \\varepsilon_0 c_0 n_0 \\operatorname{Re}(n_0) n_2 |E|^2 E + + 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. 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. Note ---- @@ -141,15 +318,103 @@ class NonlinearSusceptibility(NonlinearSpec): Example ------- - >>> medium = Medium(permittivity=2, nonlinear_spec=NonlinearSusceptibility(chi3=1)) + >>> kerr_model = KerrNonlinearity(n2=1) """ - chi3: float = pd.Field( - ..., title="Chi3", description="Chi3 nonlinear susceptibility.", units="um^2 / V^2" + n2: Complex = pd.Field( + 0, + title="Nonlinear refractive index", + description="Nonlinear refractive index in the Kerr nonlinearity.", + units=f"{MICROMETER}^2 / {WATT}", + ) + + n0: Optional[Complex] = pd.Field( + None, + title="Complex linear refractive index", + description="Complex linear refractive index of the medium, computed for instance using " + "'medium.nk_model'. If not provided, it is calculated automatically using the central " + "frequencies of the simulation sources (as long as these are all equal).", ) + 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 medium.allow_gain: + chi_imag = np.imag(n2 * n0 * np.real(n0)) + if chi_imag < 0: + raise ValidationError( + "For passive medium, 'n2' in 'KerrNonlinearity' must satisfy " + f"'Im(n2 * n0 * Re(n0)) >= 0'. Currently, this quantity equals '{chi_imag}', " + f"and the linear index is 'n0={n0}'. To simulate gain medium, please set " + "'allow_gain=True' in the medium class. Caution: simulations containing " + "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, []) + + +NonlinearModelType = Union[NonlinearSusceptibility, TwoPhotonAbsorption, KerrNonlinearity] + + +class NonlinearSpec(ABC, Tidy3dBaseModel): + """Abstract specification for adding nonlinearities to a medium. + + Note + ---- + The nonlinear constitutive relation is solved iteratively; it may not converge + for strong nonlinearities. Increasing ``num_iters`` can help with convergence. -NonlinearSpecType = Union[NonlinearSusceptibility] + Example + ------- + >>> nonlinear_susceptibility = NonlinearSusceptibility(chi3=1) + >>> nonlinear_spec = NonlinearSpec(models=[nonlinear_susceptibility]) + >>> medium = Medium(permittivity=2, nonlinear_spec=nonlinear_spec) + """ + + models: Tuple[NonlinearModelType, ...] = pd.Field( + (), + title="Nonlinear models", + description="The nonlinear models present in this nonlinear spec. " + "Nonlinear models of different types are additive. " + "Multiple nonlinear models of the same type are not allowed.", + ) + + num_iters: pd.PositiveInt = pd.Field( + NONLINEAR_DEFAULT_NUM_ITERS, + title="Number of iterations", + description="Number of iterations for solving nonlinear constitutive relation.", + ) + + @pd.validator("models", always=True) + def _no_duplicate_models(cls, val): + """Ensure each type of model appears at most once.""" + if val is None: + return val + models = [model.__class__ for model in val] + models_unique = set(models) + if len(models) != len(models_unique): + raise ValidationError( + "Multiple 'NonlinearModels' of the same type " + "were found in a single 'NonlinearSpec'. Please ensure that " + "each type of 'NonlinearModel' appears at most once in a single 'NonlinearSpec'." + ) + return val + + @pd.validator("num_iters", always=True) + def _validate_num_iters(cls, val, values): + """Check that num_iters is not too large.""" + if val > NONLINEAR_MAX_NUM_ITERS: + raise ValidationError( + "'NonlinearSpec.num_iters' must be less than " + f"{NONLINEAR_MAX_NUM_ITERS}, currently {val}." + ) + return val class AbstractMedium(ABC, Tidy3dBaseModel): @@ -174,7 +439,7 @@ class AbstractMedium(ABC, Tidy3dBaseModel): "useful in some cases.", ) - nonlinear_spec: NonlinearSpecType = pd.Field( + nonlinear_spec: Union[NonlinearSpec, NonlinearSusceptibility] = pd.Field( None, title="Nonlinear Spec", description="Nonlinear spec applied on top of the base medium properties.", @@ -186,15 +451,55 @@ class AbstractMedium(ABC, Tidy3dBaseModel): description="Modulation spec applied on top of the base medium properties.", ) - @pd.validator("nonlinear_spec", always=True) - def _validate_nonlinear_spec(cls, val): + @cached_property + def _nonlinear_models(self) -> NonlinearSpec: + """The nonlinear models in the nonlinear_spec.""" + if self.nonlinear_spec is None: + return [] + if isinstance(self.nonlinear_spec, NonlinearModel): + return [self.nonlinear_spec] + if self.nonlinear_spec.models is None: + return [] + return self.nonlinear_spec.models + + @cached_property + def _nonlinear_num_iters(self) -> pd.PositiveInt: + """The num_iters of the nonlinear_spec.""" + if self.nonlinear_spec is None: + return 0 + if isinstance(self.nonlinear_spec, NonlinearModel): + if self.nonlinear_spec.numiters is None: + return 1 # old default value for backwards compatibility + return self.nonlinear_spec.numiters + return self.nonlinear_spec.num_iters + + def _post_init_validators(self) -> None: + """Call validators taking `self` that get run after init.""" + self._validate_nonlinear_spec() + + def _validate_nonlinear_spec(self): """Check compatibility with nonlinear_spec.""" - if val is None: - return val - raise ValidationError( - f"A 'nonlinear_spec' of class {type(val)} is not " - f"currently supported for medium class {cls}." - ) + if self.nonlinear_spec is None: + return + if isinstance(self.nonlinear_spec, NonlinearModel): + log.warning( + "The API for 'nonlinear_spec' has changed. " + "The old usage 'nonlinear_spec=model' is deprecated and will be removed " + "in a future release. The new usage is " + r"'nonlinear_spec=NonlinearSpec(models=\[model])'." + ) + for model in self._nonlinear_models: + model._validate_medium_type(self) + model._validate_medium(self) + if ( + isinstance(self.nonlinear_spec, NonlinearSpec) + and isinstance(model, NonlinearSusceptibility) + and model.numiters is not None + ): + raise ValidationError( + "'NonlinearSusceptibility.numiters' is deprecated. " + "Please use 'NonlinearSpec.num_iters' instead." + ) heat_spec: Optional[HeatSpecType] = pd.Field( None, @@ -681,16 +986,6 @@ class Medium(AbstractMedium): units=CONDUCTIVITY, ) - @pd.validator("nonlinear_spec", always=True) - def _validate_nonlinear_spec(cls, val): - """Check compatibility with nonlinear_spec.""" - if val is None or isinstance(val, NonlinearSusceptibility): - return val - raise ValidationError( - f"A 'nonlinear_spec' of class {type(val)} is not " - f"currently supported for medium class {cls}." - ) - @pd.validator("conductivity", always=True) def _passivity_validation(cls, val, values): """Assert passive medium if `allow_gain` is False.""" @@ -810,16 +1105,6 @@ class CustomIsotropicMedium(AbstractCustomMedium, Medium): units=CONDUCTIVITY, ) - @pd.validator("nonlinear_spec", always=True) - def _validate_nonlinear_spec(cls, val): - """Check compatibility with nonlinear_spec.""" - if val is None: - return val - raise ValidationError( - f"A 'nonlinear_spec' of class {type(val)} is not " - f"currently supported for medium class {cls}." - ) - @pd.validator("permittivity", always=True) def _eps_inf_greater_no_less_than_one(cls, val): """Assert any eps_inf must be >=1""" diff --git a/tidy3d/components/simulation.py b/tidy3d/components/simulation.py index 8b3fbae01..ecb635187 100644 --- a/tidy3d/components/simulation.py +++ b/tidy3d/components/simulation.py @@ -26,7 +26,7 @@ from .boundary import PML, StablePML, Absorber, AbsorberSpec from .structure import Structure from .source import SourceType, PlaneWave, GaussianBeam, AstigmaticGaussianBeam, CustomFieldSource -from .source import CustomCurrentSource, CustomSourceTime +from .source import CustomCurrentSource, CustomSourceTime, ContinuousWave from .source import TFSF, Source, ModeSource from .monitor import MonitorType, Monitor, FreqMonitor, SurfaceIntegrationMonitor from .monitor import AbstractModeMonitor, FieldMonitor @@ -821,6 +821,21 @@ def _check_normalize_index(cls, val, values): if sources[val].source_time.amplitude == 0: raise ValidationError("Cannot set 'normalize_index' to source with zero amplitude.") + # Warn if normalizing by a ContinuousWave or CustomSourceTime source, if frequency-domain monitors are present. + if isinstance(sources[val].source_time, ContinuousWave): + log.warning( + f"'normalize_index' {val} is a source with 'ContinuousWave' " + "time dependence. Normalizing frequency-domain monitors by this " + "source is not meaningful because field decay does not occur. " + "Consider setting 'normalize_index' to 'None' instead." + ) + if isinstance(sources[val].source_time, CustomSourceTime): + log.warning( + f"'normalize_index' {val} is a source with 'CustomSourceTime' " + "time dependence. Normalizing frequency-domain monitors by this " + "source is only meaningful if field decay occurs." + ) + return val """ Post-init validators """ @@ -830,6 +845,7 @@ def _post_init_validators(self) -> None: _ = self.scene self._validate_no_structures_pml() self._validate_tfsf_nonuniform_grid() + self._validate_nonlinear_specs() def _validate_no_structures_pml(self) -> None: """Ensure no structures terminate / have bounds inside of PML.""" @@ -901,6 +917,15 @@ def _validate_tfsf_nonuniform_grid(self) -> None: custom_loc=["sources", source_ind], ) + def _validate_nonlinear_specs(self) -> None: + """Run :class:`.NonlinearSpec` validators that depend on knowing the central + frequencies of the sources.""" + freqs = np.array([source.source_time.freq0 for source in self.sources]) + for medium in self.scene.mediums: + if medium.nonlinear_spec is not None: + for model in medium._nonlinear_models: + model._validate_medium_freqs(medium, freqs) + """ Pre submit validation (before web.upload()) """ def validate_pre_upload(self, source_required: bool = True) -> None: diff --git a/tidy3d/components/source.py b/tidy3d/components/source.py index 8d669d8d5..ef91b1b99 100644 --- a/tidy3d/components/source.py +++ b/tidy3d/components/source.py @@ -168,6 +168,11 @@ class ContinuousWave(Pulse): """Source time dependence that ramps up to continuous oscillation and holds until end of simulation. + Note + ---- + Field decay will not occur, so the simulation will run for the full ``run_time``. + Also, source normalization of frequency-domain monitors is not meaningful. + Example ------- >>> cw = ContinuousWave(freq0=200e12, fwidth=20e12) @@ -200,6 +205,13 @@ class CustomSourceTime(Pulse): e^{i \\cdot phase - 2 \\pi i \\cdot freq0 \\cdot t} \\cdot \\ envelope(t - offset / (2 \\pi \\cdot fwidth)) + Note + ---- + Depending on the envelope, field decay may not occur. + If field decay does not occur, then the simulation will run for the full ``run_time``. + Also, if field decay does not occur, then source normalization of frequency-domain + monitors is not meaningful. + Note ---- The source time dependence is linearly interpolated to the simulation time steps. diff --git a/tidy3d/constants.py b/tidy3d/constants.py index 96429096e..aeacf056e 100644 --- a/tidy3d/constants.py +++ b/tidy3d/constants.py @@ -47,6 +47,8 @@ KELVIN = "K" CMCUBE = "cm^3" PERCMCUBE = "1/cm^3" +WATT = "W" +VOLT = "V" THERMAL_CONDUCTIVITY = "W/(um*K)" SPECIFIC_HEAT_CAPACITY = "J/(kg*K)"