Added two-photon absorption nonlinearity

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.

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,

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])

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,43 +269,51 @@ 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)
 
     # 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 * 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])
     assert np.allclose(cst.amp_time(sim.tmesh[0]), [0], rtol=0, atol=ATOL)
     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)

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]

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",