fix Medium.conductivity gradients

CHANGELOG.md

@@ -35,6 +35,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Error when plotting mode plane PML and the simulation has symmetry.
 - Validators using `TriangleMesh.intersections_plane` will fall back on bounding box in case the method fails for a non-watertight mesh.
 - Bug when running the same `ModeSolver` first locally then remotely, or vice versa, in which case the cached data from the first run is always returned.
+- Gradient inaccuracies in `PolySlab.vertices`, `Medium.conductivity`, and `DiffractionData` s-polarization.
+- Adjoint field monitors no longer store H fields, which aren't needed for gradient calculation.
+- `MeshOverrideStructures` in a `Simulation.GridSpec` are properly handled to remove any derivative tracers.
 
 ## [2.7.1] - 2024-07-10
 

tests/test_components/test_autograd.py

@@ -36,7 +36,7 @@
 TEST_POLYSLAB_SPEED = False
 
 
-TEST_MODES = ("pipeline", "adjoint", "numerical", "speed")
+TEST_MODES = ("pipeline", "adjoint", "speed")
 TEST_MODE = "speed" if TEST_POLYSLAB_SPEED else "pipeline"
 
 # number of elements in the parameters / input to the objective function
@@ -48,7 +48,7 @@
 params0 /= np.linalg.norm(params0)
 
 # whether to plot the simulation within the objective function
-PLOT_SIM = True
+PLOT_SIM = False
 
 # whether to include a call to `objective(params)` in addition to gradient
 CALL_OBJECTIVE = False
@@ -74,7 +74,7 @@
 
 # shape of the custom medium
 DA_SHAPE_X = 1 if IS_3D else 1
-DA_SHAPE = (DA_SHAPE_X, 1_0, 1_0) if TEST_CUSTOM_MEDIUM_SPEED else (DA_SHAPE_X, 12, 12)
+DA_SHAPE = (DA_SHAPE_X, 1_000, 1_000) if TEST_CUSTOM_MEDIUM_SPEED else (DA_SHAPE_X, 12, 12)
 
 # number of vertices in the polyslab
 NUM_VERTICES = 100_000 if TEST_POLYSLAB_SPEED else 30
@@ -105,13 +105,13 @@
         )
     ],
     structures=[
-        # td.Structure(
-        #     geometry=td.Box(
-        #         size=(0.5, 0.5, LZ / 2),
-        #         center=(0, 0, LZ / 2),
-        #     ),
-        #     medium=td.Medium(permittivity=1.0),
-        # )
+        td.Structure(
+            geometry=td.Box(
+                size=(0.5, 0.5, LZ / 2),
+                center=(0, 0, LZ / 2),
+            ),
+            medium=td.Medium(permittivity=2.0),
+        )
     ],
     monitors=[
         td.FieldMonitor(
@@ -254,11 +254,15 @@ def make_structures(params: anp.ndarray) -> dict[str, td.Structure]:
 
     # Structure with variable .medium
     eps = 1 + anp.abs(vector @ params)
+    sigma = 0.1 * (anp.tanh(vector @ params) + 1)
+
+    permittivity, conductivity = eps, sigma
+    # permittivity, conductivity = eps, 0.0 # only permittivity
+    # permittivity, conductivity = 2.0, sigma # only conductivity
 
-    conductivity = 0  # eps / 10.0
     medium = td.Structure(
         geometry=box,
-        medium=td.Medium(permittivity=eps, conductivity=conductivity),
+        medium=td.Medium(permittivity=permittivity, conductivity=conductivity),
     )
 
     # Structure with variable Box.center
@@ -336,8 +340,8 @@ def make_structures(params: anp.ndarray) -> dict[str, td.Structure]:
             vertices=vertices,
             slab_bounds=(-0.5, 0.5),
             axis=0,
-            sidewall_angle=0.0,  # 1,
-            dilation=0.00,  # 1,
+            sidewall_angle=0.01,
+            dilation=0.01,
         ),
         medium=med,
     )
@@ -424,16 +428,6 @@ def mode_postprocess_fn(sim_data, mnt_data):
     )
 
     def diff_postprocess_fn(sim_data, mnt_data):
-        # with minus sign on amplitude
-        # orders_y= 0  (RMS 0.0615)
-        # orders_y= 1  (RMS 0.0177)
-        # orders_y= -1  (RMS 0.02)
-        # orders_y= 2  (RMS 2)
-        # orders_y= 3  (RMS 0.3965)
-        # all (RMS 1.965)
-
-        # import pdb; pdb.set_trace()
-
         return anp.sum(abs(mnt_data.amps.sel(polarization=["s", "p"]).values) ** 2)
 
     field_vol = td.FieldMonitor(
@@ -448,10 +442,10 @@ def field_vol_postprocess_fn(sim_data, mnt_data):
         for _, val in mnt_data.field_components.items():
             value = value + abs(anp.sum(val.values))
         # field components numerical is 3x higher
-        # intensity = anp.nan_to_num(anp.sum(sim_data.get_intensity(mnt_data.monitor.name).values))
-        # value += intensity
+        intensity = anp.nan_to_num(anp.sum(sim_data.get_intensity(mnt_data.monitor.name).values))
+        value += intensity
         # intensity numerical is 4.79x higher
-        # value += anp.sum(mnt_data.flux.values)
+        value += anp.sum(mnt_data.flux.values)
         # flux is 18.4x lower
         return value
 
@@ -464,8 +458,8 @@ def field_vol_postprocess_fn(sim_data, mnt_data):
 
     def field_point_postprocess_fn(sim_data, mnt_data):
         value = 0.0
-        # for _, val in mnt_data.field_components.items():
-        #     value += abs(anp.sum(val.values))
+        for _, val in mnt_data.field_components.items():
+            value += abs(anp.sum(val.values))
         value += anp.sum(sim_data.get_intensity(mnt_data.monitor.name).values)
         return value
 
@@ -591,15 +585,12 @@ def test_polyslab_axis_ops(axis):
     basis_vecs = p.edge_basis_vectors(edges=edges)
 
 
-@pytest.mark.parametrize("structure_key, monitor_key", (("medium", "field_vol"),))
-def test_autograd_numerical(structure_key, monitor_key):
+@pytest.mark.parametrize("structure_key, monitor_key", (("medium", "mode"),))
+def _test_autograd_numerical(structure_key, monitor_key):
     """Test an objective function through tidy3d autograd."""
 
     import tidy3d.web as web
 
-    if False and TEST_MODE != "numerical":
-        return
-
     fn_dict = get_functions(structure_key, monitor_key)
     make_sim = fn_dict["sim"]
     postprocess = fn_dict["postprocess"]
@@ -618,7 +609,7 @@ def objective(*args):
     assert anp.all(grad != 0.0), "some gradients are 0"
 
     # numerical gradients
-    delta = 4e-3
+    delta = 1e-1
     sims_numerical = {}
 
     params_num = np.zeros((N_PARAMS, N_PARAMS))
@@ -666,10 +657,6 @@ def task_name_fn(i: int, sign: int) -> str:
     print(f"|grad| / |grad_num| = {norm_factor:.4f}")
     print(f"avg(diff(objectives)) = {diff_objectives_num:.4f}")
 
-    # import pdb
-
-    # pdb.set_trace()
-
 
 @pytest.mark.parametrize("structure_key, monitor_key", (("polyslab", "mode"),))
 def test_autograd_objective_tyler(use_emulated_run, structure_key, monitor_key):
@@ -1228,7 +1215,7 @@ def get_amps(sim_data: td.SimulationData, mnt_name: str) -> xr.DataArray:
 
 def power(amps: xr.DataArray) -> float:
     """Reduce a selected DataArray into just a float for objective function."""
-    return anp.sum(anp.abs(amps.sel(orders_x=[0], orders_y=[0]) ** 2))
+    return anp.sum(anp.abs(amps.values) ** 2)
 
 
 def postprocess_0_src(sim_data: td.SimulationData) -> float:

tidy3d/components/data/monitor_data.py

@@ -1098,7 +1098,7 @@ def shift_value(coords) -> float:
             for name, field_component in self.field_components.items():
                 field_component = field_component.sel(f=freq0)
                 forward_amps = field_component.values
-                values = -1j * 1 / 3.0 * forward_amps
+                values = -1j * forward_amps
                 coords = dict(field_component.coords.copy())
                 for dim, key in enumerate("xyz"):
                     coords[key] = np.array(coords[key]) - source_geo.center[dim]

tidy3d/components/medium.py

@@ -1636,7 +1636,10 @@ def derivative_eps_sigma_volume(
         freqs = vjp_eps_complex.coords["f"].values
         values = vjp_eps_complex.values
 
-        eps_vjp, sigma_vjp = self.eps_complex_to_eps_sigma(eps_complex=values, freq=freqs)
+        # vjp of eps_complex_to_eps_sigma
+        omegas = 2 * np.pi * freqs
+        eps_vjp = np.real(values)
+        sigma_vjp = -np.imag(values) / omegas / EPSILON_0
 
         eps_vjp = np.sum(eps_vjp)
         sigma_vjp = np.sum(sigma_vjp)

tidy3d/web/api/autograd/autograd.py

@@ -595,22 +595,8 @@ def vjp(data_fields_vjp: AutogradFieldMap) -> AutogradFieldMap:
         task_name_adj = str(task_name) + "_adjoint"
 
         if local_gradient:
-            sim_adj = sim_adj.updated_copy(
-                monitors=list(sim_adj.monitors)
-                + [
-                    td.FieldMonitor(
-                        size=(td.inf, td.inf, td.inf),
-                        center=(0, 0, 0),
-                        freqs=[299792458000000.0],
-                        name="adjoint fields",
-                    )
-                ]
-            )
-
             sim_data_adj, _ = _run_tidy3d(sim_adj, task_name=task_name_adj, **run_kwargs)
 
-            # import pdb; pdb.set_trace()
-
             vjp_traced_fields = postprocess_adj(
                 sim_data_adj=sim_data_adj,
                 sim_data_orig=sim_data_orig,