[FXC-5526] add non-colocated dot and outer_dot for improved accuracy

[dmarek-flex]

Note

Medium Risk
Touches core numerical integration paths for overlaps/flux and changes mode normalization, which can subtly affect results and downstream computations; mitigated by added coverage for Yee-area consistency, non-uniform grids, and degenerate-mode orthogonality.

Overview
Improves FieldData/ModeData overlap calculations by adding a non-colocated path for dot/outer_dot that avoids boundary colocation/interpolation and instead integrates using Yee-staggered differential area elements (cell×dual and dual×cell) for better accuracy on non-uniform grids.

Adds NumPy utilities (_dot_numpy, _outer_dot_numpy, _complex_power_flow_numpy, _instantaneous_power_flow_numpy, plus coordinate broadcast/intersection helpers) and refactors dot/outer_dot to use these fast paths, returning FreqDataArray when no mode_index is present and adding a use_colocated_fields override.

Extends flux handling for non-colocated monitors via _complex_flux_non_colocated (freq-domain) and _flux_non_colocated (time-domain), wires this into the microwave impedance calculator’s fallback-to-flux logic, and updates mode normalization to use self-overlap (self.dot(self)) rather than sqrt(|flux|).

Updates the mode solver’s internal overlap routines to reuse the new NumPy overlap kernels, and adds/adjusts tests to validate Yee area element construction (including 2D unit handling and non-uniform meshes) and to make mode-solver colocation comparisons normalization-invariant while tightening degenerate-mode orthogonality tolerances.

^{Written by Cursor Bugbot for commit 15dc89d. This will update automatically on new commits. Configure here.}

[cursor]

Cursor Bugbot has reviewed your changes and found 1 potential issue.

^{Bugbot Autofix is OFF. To automatically fix reported issues with Cloud Agents, enable Autofix in the Cursor dashboard.}

[cursor]

Paired numpy advanced indexing corrupts array selection

Medium Severity

When both f_sel and m_sel in the indexer are numpy arrays (not slices), arr[tuple(indexer)] performs paired advanced indexing instead of independent (outer) indexing. This happens when the other dataset has size-1 frequency AND size-1 mode_index while self has different counts of each (N_f ≠ N_m), causing either an IndexError from incompatible broadcast shapes, or silently wrong results (paired selection instead of outer product) when N_f == N_m. _get_broadcast_selection returns np.zeros(N_f) and np.zeros(N_m) arrays for the two dimensions, and placing both in the same indexer tuple triggers numpy's coupled advanced indexing semantics. The fix likely requires applying selections sequentially or using np.ix_.

Additional Locations (1)

• tidy3d/components/data/utils.py#L122-L129

 

[github-actions]

Diff Coverage

Diff: origin/develop...HEAD, staged and unstaged changes

• tidy3d/components/data/monitor_data.py (93.2%): Missing lines 465,533,576,583,811,891-894,901,947,966,1154,1160,1175,1399,1405,1420,2052,2062,2078
• tidy3d/components/data/utils.py (91.9%): Missing lines 135,137-138,223-224,226
• tidy3d/components/microwave/impedance_calculator.py (92.3%): Missing lines 141
• tidy3d/components/mode/solver.py (100%)

Summary

• Total: 405 lines
• Missing: 28 lines
• Coverage: 93%

tidy3d/components/data/monitor_data.py

Lines 461-469

  461     def _normal_dim(self) -> str:
  462         """For a 2D monitor data, return the name of the normal dimension. Raise if cannot
  463         confirm that the associated monitor is 2D."""
  464         if len(self.monitor.zero_dims) != 1:
! 465             raise DataError("Data must be 2D to get normal dimension.")
  466         normal_dim = "xyz"[self.monitor.size.index(0)]
  467         return normal_dim
  468 
  469     @property

---

Lines 529-537

  529             area elements for Yee grid integration. Each DataArray has dimensions
  530             corresponding to the two tangential dimensions of the monitor plane.
  531         """
  532         if not self.grid_expanded:
! 533             raise DataError(
  534                 "Monitor data requires 'grid_expanded' to compute Yee grid integration sizes."
  535             )
  536 
  537         _, plane_inds = self.monitor.pop_axis([0, 1, 2], self.monitor.size.index(0.0))

---

Lines 572-580

  572             # Determine integration bounds
  573             if truncate_to_monitor_bounds:
  574                 # Use monitor bounds, handling inf by finding grid boundary outside field data
  575                 if np.isinf(mnt_min):
! 576                     integration_min = self._find_enclosing_boundary(
  577                         e_field_centers[0], full_boundaries, "lower"
  578                     )
  579                 else:
  580                     integration_min = mnt_min

---

Lines 579-587

  579                 else:
  580                     integration_min = mnt_min
  581 
  582                 if np.isinf(mnt_max):
! 583                     integration_max = self._find_enclosing_boundary(
  584                         e_field_centers[-1], full_boundaries, "upper"
  585                     )
  586                 else:
  587                     integration_max = mnt_max

---

Lines 807-815

  807         field_data = self.symmetry_expanded_copy
  808         if len(self.monitor.zero_dims) != 1:
  809             return field_data
  810 
! 811         normal_dim = self._normal_dim
  812         update = {"grid_primal_correction": 1.0, "grid_dual_correction": 1.0}
  813         for field_name, field in field_data.field_components.items():
  814             eig_val = self.symmetry_eigenvalues[field_name](normal_dim)
  815             if eig_val < 0:

---

Lines 887-898

  887             if i != u_axis and i != v_axis and (normal_axis is None or i != normal_axis)
  888         ]
  889         # Adjust axes for removed normal dimension
  890         if normal_axis is not None:
! 891             adjusted_other = [a if a < normal_axis else a - 1 for a in other_axes]
! 892             adjusted_u = u_axis if u_axis < normal_axis else u_axis - 1
! 893             adjusted_v = v_axis if v_axis < normal_axis else v_axis - 1
! 894             transpose_axes = [*adjusted_other, adjusted_u, adjusted_v]
  895         else:
  896             transpose_axes = [*other_axes, u_axis, v_axis]
  897 
  898         def prepare_numpy(arr: np.ndarray) -> np.ndarray:

---

Lines 897-905

  897 
  898         def prepare_numpy(arr: np.ndarray) -> np.ndarray:
  899             """Prepare array: squeeze normal dim if present, then transpose."""
  900             if normal_axis is not None:
! 901                 arr = np.squeeze(arr, axis=normal_axis)
  902             return arr.transpose(transpose_axes)
  903 
  904         prepped_fields = {key: prepare_numpy(field.values) for key, field in fields.items()}

---

Lines 943-951

  943         falls back to the standard ``complex_flux``.
  944         """
  945 
  946         if self.monitor.colocate:
! 947             raise ValueError(
  948                 "Non-colocated flux is only available for field data from monitors with 'colocate=False'."
  949             )
  950 
  951         fields = self._tangential_fields

---

Lines 962-970

  962         flux_result = _complex_power_flow_numpy(E, H, dS_numpy)
  963 
  964         if "mode_index" in final_coords:
  965             return FreqModeDataArray(flux_result, coords=final_coords)
! 966         return FluxDataArray(flux_result, coords=final_coords)
  967 
  968     @cached_property
  969     def mode_area(self) -> FreqModeDataArray:
  970         r"""Effective mode area corresponding to a 2D monitor.

---

Lines 1150-1158

  1150             arr = arr[tuple(indexer)]
  1151 
  1152             # Squeeze out normal dimension if present
  1153             if normal_axis is not None:
! 1154                 arr = np.squeeze(arr, axis=normal_axis)
  1155 
  1156             # Compute working axes after potential normal squeeze
  1157             def adjust_axis(ax: Optional[int]) -> Optional[int]:
  1158                 if normal_axis is None or ax is None:

---

Lines 1156-1164

  1156             # Compute working axes after potential normal squeeze
  1157             def adjust_axis(ax: Optional[int]) -> Optional[int]:
  1158                 if normal_axis is None or ax is None:
  1159                     return ax
! 1160                 return ax if ax < normal_axis else ax - 1
  1161 
  1162             f_ax = adjust_axis(f_axis)
  1163             m_ax = adjust_axis(mode_axis) if has_mode else None
  1164             u_ax = adjust_axis(u_axis)

---

Lines 1171-1179

  1171                 m_ax = 1
  1172                 if f_ax >= 1:
  1173                     f_ax += 1
  1174                 if u_ax >= 1:
! 1175                     u_ax += 1
  1176                 if v_ax >= 1:
  1177                     v_ax += 1
  1178 
  1179             # Transpose to (f, mode_index, u, v)

---

Lines 1395-1403

  1395             arr = arr[tuple(indexer)]
  1396 
  1397             # Squeeze out normal dimension if present
  1398             if normal_axis is not None:
! 1399                 arr = np.squeeze(arr, axis=normal_axis)
  1400 
  1401             # Compute working axes after potential normal squeeze
  1402             def adjust_axis(ax: Optional[int]) -> Optional[int]:
  1403                 if normal_axis is None or ax is None:

---

Lines 1401-1409

  1401             # Compute working axes after potential normal squeeze
  1402             def adjust_axis(ax: Optional[int]) -> Optional[int]:
  1403                 if normal_axis is None or ax is None:
  1404                     return ax
! 1405                 return ax if ax < normal_axis else ax - 1
  1406 
  1407             f_ax = adjust_axis(f_axis)
  1408             m_ax = adjust_axis(mode_axis) if has_mode else None
  1409             u_ax = adjust_axis(u_axis)

---

Lines 1416-1424

  1416                 m_ax = 1
  1417                 if f_ax >= 1:
  1418                     f_ax += 1
  1419                 if u_ax >= 1:
! 1420                     u_ax += 1
  1421                 if v_ax >= 1:
  1422                     v_ax += 1
  1423 
  1424             # Transpose to (f, mode_index, u, v)

---

Lines 2048-2056

  2048 
  2049     @cached_property
  2050     def complex_flux(self) -> DataArray:
  2051         """Complex flux is not defined for time-domain data."""
! 2052         raise DataError("Complex power flow is not defined for time-domain data.")
  2053 
  2054     def _flux_non_colocated(self) -> FluxTimeDataArray:
  2055         """Non-colocated complex flux for internal use.

---

Lines 2058-2066

  2058         flux integration when ``self.monitor.colocate is False``.
  2059         """
  2060 
  2061         if self.monitor.colocate:
! 2062             raise ValueError(
  2063                 "Non-colocated flux is only available for field data from monitors with 'colocate=False'."
  2064             )
  2065         dim1, dim2 = self._tangential_dims
  2066         tangential_dims = self._tangential_dims

---

Lines 2074-2082

  2074         # Prepare arrays: remove normal dim if present, put spatial dims last
  2075         def prepare(data_array: DataArray) -> DataArray:
  2076             prepared = data_array
  2077             if normal_dim in prepared.dims:
! 2078                 prepared = prepared.squeeze(dim=normal_dim, drop=True)
  2079             return prepared.transpose(..., *tangential_dims)
  2080 
  2081         # Extract and prepare field components
  2082         Eu = np.real(prepare(fields["E" + dim1]).to_numpy())

---

tidy3d/components/data/utils.py

Lines 131-142

  131         # Fast path: exact match, use slices (creates views, not copies)
  132         return slice(None), slice(None), vals_self
  133     else:
  134         # Intersection
! 135         common, idx_self, idx_other = np.intersect1d(vals_self, vals_other, return_indices=True)
  136         # Preserve order from self
! 137         order = np.argsort(idx_self)
! 138         return idx_self[order], idx_other[order], common[order]
  139 
  140 
  141 def _get_intersection_selection(
  142     vals_self: np.ndarray, vals_other: np.ndarray

---

Lines 219-230

  219     else:
  220         # Non-bidirectional: 0.5 * integral(E1 x H2) dS
  221         # At Eu/Hv location: E1u * H2v
  222         # At Ev/Hu location: E1v * H2u
! 223         term_EuHv = (E1u * H2v) * dS_EuHv
! 224         term_EvHu = (E1v * H2u) * dS_EvHu
  225         # Sum over spatial dimensions (last two)
! 226         return 0.5 * np.sum(term_EuHv - term_EvHu, axis=(-2, -1))
  227 
  228 
  229 def _outer_dot_numpy(
  230     E1: tuple[np.ndarray, np.ndarray],

---

tidy3d/components/microwave/impedance_calculator.py

Lines 137-145

  137             if isinstance(em_field.monitor, ModeMonitor):
  138                 flux_sign = 1 if em_field.monitor.store_fields_direction == "+" else -1
  139             if isinstance(em_field, FieldTimeData):
  140                 if em_field.monitor.colocate:
! 141                     flux = flux_sign * em_field.flux
  142                 else:
  143                     flux = flux_sign * em_field._flux_non_colocated()
  144             else:
  145                 if em_field.monitor.colocate:

---