From 7a5ea1c218ae580f86dff6e824f81695d558588a Mon Sep 17 00:00:00 2001 From: Momchil Minkov Date: Mon, 15 Dec 2025 10:03:01 +0100 Subject: [PATCH 1/3] feat: charge solver adaptive mesh based on doping --- tidy3d/__init__.py | 2 + tidy3d/components/material/tcad/charge.py | 165 +++++++++++++++++- tidy3d/components/tcad/doping.py | 80 +++++++++ tidy3d/components/tcad/grid.py | 80 +++++++++ .../components/tcad/simulation/heat_charge.py | 158 +++++++++++++++++ 5 files changed, 484 insertions(+), 1 deletion(-) diff --git a/tidy3d/__init__.py b/tidy3d/__init__.py index 18f7d99047..11b8de9f4b 100644 --- a/tidy3d/__init__.py +++ b/tidy3d/__init__.py @@ -92,6 +92,7 @@ from tidy3d.components.tcad.generation_recombination import FossumCarrierLifetime from tidy3d.components.tcad.grid import ( DistanceUnstructuredGrid, + DopingGradientRefinementSpec, GridRefinementLine, GridRefinementRegion, UniformUnstructuredGrid, @@ -594,6 +595,7 @@ def set_logging_level(level: str) -> None: "DirectivityMonitor", "DistanceUnstructuredGrid", "DistributedGeneration", + "DopingGradientRefinementSpec", "Drude", "DualValleyEffectiveDOS", "EMECoefficientData", diff --git a/tidy3d/components/material/tcad/charge.py b/tidy3d/components/material/tcad/charge.py index 94d8cdb778..fec0e0bf8f 100644 --- a/tidy3d/components/material/tcad/charge.py +++ b/tidy3d/components/material/tcad/charge.py @@ -2,8 +2,9 @@ from __future__ import annotations -from typing import Union +from typing import TYPE_CHECKING, Union +import numpy as np import pydantic.v1 as pd from tidy3d.components.data.data_array import SpatialDataArray @@ -21,6 +22,10 @@ from tidy3d.constants import CONDUCTIVITY, ELECTRON_VOLT, PERCMCUBE, PERMITTIVITY from tidy3d.log import log +if TYPE_CHECKING: + from tidy3d.components.tcad.grid import DopingGradientRefinementSpec, GridRefinementRegion + from tidy3d.components.types import Bound + class AbstractChargeMedium(AbstractMedium): """Abstract class for Charge specifications @@ -375,3 +380,161 @@ def check_na_uses_model(cls, val, values): ) return (ConstantDoping(concentration=val),) return val + + def compute_doping_refinement_regions( + self, + bounds: Bound, + spec: DopingGradientRefinementSpec, + normal_axis: int, + ) -> list[GridRefinementRegion]: + """Compute mesh refinement regions based on doping gradients. + + This method samples the net doping (|N_d - N_a|) on a 2D grid within + the given bounds, computes the gradient magnitude, and generates + :class:`GridRefinementRegion` objects in areas of rapid doping variation. + + Parameters + ---------- + bounds : Bound + Tuple of (min_corner, max_corner) defining the sampling region. + spec : DopingGradientRefinementSpec + Specification controlling refinement parameters. + normal_axis : int + The axis normal to the 2D plane (0=x, 1=y, 2=z). + + Returns + ------- + list[GridRefinementRegion] + List of refinement regions to be added to the mesh specification. + """ + from tidy3d.components.tcad.grid import GridRefinementRegion + + refinement_regions = [] + + # Determine the in-plane axes + axes = [0, 1, 2] + axes.remove(normal_axis) + ax1, ax2 = axes + + # Extract bounds + min_corner, max_corner = bounds + normal_pos = (min_corner[normal_axis] + max_corner[normal_axis]) / 2 + + # Create sampling grid + num_samples = spec.num_samples + coord_names = ["x", "y", "z"] + + # Create 1D coordinate arrays for the 2D plane + coords_1 = np.linspace(min_corner[ax1], max_corner[ax1], num_samples) + coords_2 = np.linspace(min_corner[ax2], max_corner[ax2], num_samples) + + # Build coords dict for sampling + coords = {coord_names[normal_axis]: np.array([normal_pos])} + coords[coord_names[ax1]] = coords_1 + coords[coord_names[ax2]] = coords_2 + + # Compute net doping and gradient for N_d + net_doping = np.zeros((num_samples, num_samples, 1)) + grad_magnitude = np.zeros((num_samples, num_samples, 1)) + + # Process donor doping (N_d) + if isinstance(self.N_d, tuple): + for doping_box in self.N_d: + doping_contrib = doping_box._get_contrib(coords, meshgrid=True) + grad_contrib = doping_box.gradient_magnitude(coords, meshgrid=True) + if doping_contrib.ndim == 2: + doping_contrib = doping_contrib[:, :, np.newaxis] + grad_contrib = grad_contrib[:, :, np.newaxis] + net_doping = net_doping + doping_contrib + grad_magnitude = np.maximum(grad_magnitude, grad_contrib) + + # Process acceptor doping (N_a) - subtract from net doping + if isinstance(self.N_a, tuple): + for doping_box in self.N_a: + doping_contrib = doping_box._get_contrib(coords, meshgrid=True) + grad_contrib = doping_box.gradient_magnitude(coords, meshgrid=True) + if doping_contrib.ndim == 2: + doping_contrib = doping_contrib[:, :, np.newaxis] + grad_contrib = grad_contrib[:, :, np.newaxis] + net_doping = net_doping - doping_contrib + grad_magnitude = np.maximum(grad_magnitude, grad_contrib) + + # Squeeze to 2D + net_doping = np.abs(net_doping.squeeze()) + grad_magnitude = grad_magnitude.squeeze() + + # Compute local mesh size: dl = N / (|grad_N| * k) + # Avoid division by zero + eps = 1e-10 + local_dl = np.where( + grad_magnitude > eps, + net_doping / (grad_magnitude * spec.resolution_factor), + spec.max_dl, + ) + + # Clamp to [min_dl, max_dl] + local_dl = np.clip(local_dl, spec.min_dl, spec.max_dl) + + # Find regions where mesh needs to be finer than max_dl + # Group contiguous fine regions into refinement boxes + needs_refinement = local_dl < spec.max_dl * 0.9 # 90% threshold + + if not np.any(needs_refinement): + return refinement_regions + + # Simple approach: create refinement regions around each connected component + # For efficiency, we use a grid-based clustering approach + from scipy import ndimage + + labeled, num_features = ndimage.label(needs_refinement) + + for label_id in range(1, num_features + 1): + mask = labeled == label_id + + # Find bounding box of this region + rows, cols = np.where(mask) + if len(rows) == 0: + continue + + row_min, row_max = rows.min(), rows.max() + col_min, col_max = cols.min(), cols.max() + + # Get the minimum dl in this region + dl_internal = local_dl[mask].min() + + # Compute region bounds in physical coordinates + region_min_1 = coords_1[row_min] + region_max_1 = coords_1[row_max] + region_min_2 = coords_2[col_min] + region_max_2 = coords_2[col_max] + + # Build center and size for the refinement region + center = [0.0, 0.0, 0.0] + size = [0.0, 0.0, 0.0] + + center[normal_axis] = normal_pos + size[normal_axis] = 0.0 # 2D refinement + + center[ax1] = (region_min_1 + region_max_1) / 2 + size[ax1] = region_max_1 - region_min_1 + + center[ax2] = (region_min_2 + region_max_2) / 2 + size[ax2] = region_max_2 - region_min_2 + + # Add some padding + padding = dl_internal * 2 + size[ax1] = size[ax1] + 2 * padding + size[ax2] = size[ax2] + 2 * padding + + transition_thickness = dl_internal * spec.transition_factor + + refinement_regions.append( + GridRefinementRegion( + center=tuple(center), + size=tuple(size), + dl_internal=dl_internal, + transition_thickness=transition_thickness, + ) + ) + + return refinement_regions diff --git a/tidy3d/components/tcad/doping.py b/tidy3d/components/tcad/doping.py index 1d82639891..439ba54f30 100644 --- a/tidy3d/components/tcad/doping.py +++ b/tidy3d/components/tcad/doping.py @@ -66,6 +66,24 @@ def _post_init_validators(self) -> None: "The doping box must be 3D. If you want a 2D doping box, please set one of the dimensions to a large or infinite size." ) + def gradient_magnitude(self, coords: dict, meshgrid: bool = True) -> np.ndarray: + """Compute the magnitude of the doping gradient at specified coordinates. + + Parameters + ---------- + coords : dict + Dictionary with keys 'x', 'y', 'z' containing coordinate arrays. + meshgrid : bool = True + If True, coordinates are treated as defining a meshgrid. + If False, coordinates are treated as parallel 1D arrays. + + Returns + ------- + np.ndarray + Magnitude of the doping gradient |∇N| at each coordinate. + """ + raise NotImplementedError("Subclasses must implement gradient_magnitude") + class ConstantDoping(AbstractDopingBox): """ @@ -100,6 +118,11 @@ def _get_contrib(self, coords: dict, meshgrid: bool = True): return contrib.squeeze() + def gradient_magnitude(self, coords: dict, meshgrid: bool = True) -> np.ndarray: + """Gradient magnitude for constant doping is zero everywhere.""" + indices_in_box, X, _, _ = self._get_indices_in_box(coords=coords, meshgrid=meshgrid) + return np.zeros(X.shape).squeeze() + class GaussianDoping(AbstractDopingBox): """Sets a gaussian doping in the specified box. For translationally invariant behavior in one dimension, the box must have infinite size in the @@ -272,6 +295,39 @@ def _get_contrib(self, coords: dict, meshgrid: bool = True): return total_contrib.squeeze() + def gradient_magnitude(self, coords: dict, meshgrid: bool = True) -> np.ndarray: + """Compute gradient magnitude for Gaussian doping using finite differences. + + For Gaussian doping, the gradient is computed numerically since the analytical + form involves products of Gaussians in each dimension. + """ + # Use finite differences on the doping profile + doping = self._get_contrib(coords, meshgrid=meshgrid) + + if meshgrid: + # Compute gradient using finite differences + dx = coords["x"][1] - coords["x"][0] if len(coords["x"]) > 1 else 1.0 + dy = coords["y"][1] - coords["y"][0] if len(coords["y"]) > 1 else 1.0 + dz = coords["z"][1] - coords["z"][0] if len(coords["z"]) > 1 else 1.0 + + # Use np.gradient for interior points with edge handling + grad_x = ( + np.gradient(doping, dx, axis=0) if doping.shape[0] > 1 else np.zeros_like(doping) + ) + grad_y = ( + np.gradient(doping, dy, axis=1) if doping.shape[1] > 1 else np.zeros_like(doping) + ) + grad_z = ( + np.gradient(doping, dz, axis=2) if doping.shape[2] > 1 else np.zeros_like(doping) + ) + + grad_mag = np.sqrt(grad_x**2 + grad_y**2 + grad_z**2) + return grad_mag.squeeze() + else: + # For non-meshgrid, we cannot compute gradient directly + # Return zeros as a fallback (gradients require structured grid) + return np.zeros_like(doping) + class CustomDoping(AbstractDopingBox): """Sets a custom doping in the specified box. @@ -338,5 +394,29 @@ def _get_contrib(self, coords: dict, meshgrid: bool = True): return contrib.squeeze() + def gradient_magnitude(self, coords: dict, meshgrid: bool = True) -> np.ndarray: + """Compute gradient magnitude for custom doping using finite differences.""" + doping = self._get_contrib(coords, meshgrid=meshgrid) + + if meshgrid: + dx = coords["x"][1] - coords["x"][0] if len(coords["x"]) > 1 else 1.0 + dy = coords["y"][1] - coords["y"][0] if len(coords["y"]) > 1 else 1.0 + dz = coords["z"][1] - coords["z"][0] if len(coords["z"]) > 1 else 1.0 + + grad_x = ( + np.gradient(doping, dx, axis=0) if doping.shape[0] > 1 else np.zeros_like(doping) + ) + grad_y = ( + np.gradient(doping, dy, axis=1) if doping.shape[1] > 1 else np.zeros_like(doping) + ) + grad_z = ( + np.gradient(doping, dz, axis=2) if doping.shape[2] > 1 else np.zeros_like(doping) + ) + + grad_mag = np.sqrt(grad_x**2 + grad_y**2 + grad_z**2) + return grad_mag.squeeze() + else: + return np.zeros_like(doping) + DopingBoxType = Union[ConstantDoping, GaussianDoping, CustomDoping] diff --git a/tidy3d/components/tcad/grid.py b/tidy3d/components/tcad/grid.py index fcf0cf6f27..8a94d01e89 100644 --- a/tidy3d/components/tcad/grid.py +++ b/tidy3d/components/tcad/grid.py @@ -83,6 +83,77 @@ class GridRefinementRegion(Box): ) +class DopingGradientRefinementSpec(Tidy3dBaseModel): + """Specification for automatic mesh refinement based on doping gradients. + + This generates :class:`GridRefinementRegion` objects in areas where doping concentration + varies rapidly, ensuring adequate mesh resolution to capture carrier distributions accurately. + + Notes + ----- + The local mesh size is computed as: + + .. math:: + + dl = \\frac{N}{|\\nabla N| \\cdot k} + + where :math:`N` is the local doping concentration, :math:`|\\nabla N|` is the gradient + magnitude, and :math:`k` is the ``resolution_factor``. The result is clamped to + ``[min_dl, max_dl]``. + + Currently only 2D simulations are supported (one dimension with zero size). + + Example + ------- + >>> spec = DopingGradientRefinementSpec(resolution_factor=5, min_dl=0.001, max_dl=0.1) + """ + + resolution_factor: pd.PositiveFloat = pd.Field( + 5.0, + title="Resolution Factor", + description="Number of mesh points desired per decade of doping variation. " + "Higher values produce finer meshes in regions of rapid doping change.", + ) + + min_dl: pd.PositiveFloat = pd.Field( + 0.001, + title="Minimum Mesh Size", + description="Absolute minimum allowed mesh size from doping-based refinement.", + units=MICROMETER, + ) + + max_dl: pd.PositiveFloat = pd.Field( + 1.0, + title="Maximum Mesh Size", + description="Maximum mesh size in regions of slow doping variation. " + "The actual mesh size may still be smaller due to other constraints.", + units=MICROMETER, + ) + + num_samples: pd.PositiveInt = pd.Field( + 50, + title="Sampling Resolution", + description="Number of sampling points per dimension when computing doping gradients. " + "Higher values improve accuracy but increase preprocessing time.", + ) + + transition_factor: pd.PositiveFloat = pd.Field( + 10.0, + title="Transition Factor", + description="Factor multiplied by ``dl_internal`` to determine ``transition_thickness`` " + "of generated refinement regions.", + ) + + @pd.validator("max_dl", always=True) + @skip_if_fields_missing(["min_dl"]) + def _max_dl_gte_min_dl(cls, val, values): + """Error if max_dl is less than min_dl.""" + min_dl = values.get("min_dl") + if val < min_dl: + raise ValidationError("'max_dl' cannot be smaller than 'min_dl'.") + return val + + class GridRefinementLine(Tidy3dBaseModel, ABC): """Refinement line for the unstructured mesh. The cell size depends on the distance from the line.""" @@ -223,6 +294,15 @@ class DistanceUnstructuredGrid(UnstructuredGrid): ) ) + auto_doping_refinement: DopingGradientRefinementSpec = pd.Field( + None, + title="Automatic Doping-Based Refinement", + description="Specification for automatic mesh refinement based on doping gradients. " + "When provided, refinement regions are automatically generated in areas of rapid " + "doping variation within semiconductor structures. Set to ``None`` to disable. " + "By default, this is enabled with sensible defaults for 2D charge simulations.", + ) + @pd.validator("distance_bulk", always=True) @skip_if_fields_missing(["distance_interface"]) def names_exist_bcs(cls, val, values): diff --git a/tidy3d/components/tcad/simulation/heat_charge.py b/tidy3d/components/tcad/simulation/heat_charge.py index 71a09071f5..a89be84639 100644 --- a/tidy3d/components/tcad/simulation/heat_charge.py +++ b/tidy3d/components/tcad/simulation/heat_charge.py @@ -53,6 +53,7 @@ ) from tidy3d.components.tcad.grid import ( DistanceUnstructuredGrid, + GridRefinementRegion, UniformUnstructuredGrid, UnstructuredGridType, ) @@ -2006,3 +2007,160 @@ def _get_ssac_frequency_and_amplitude(self): if isinstance(bc.condition.source, SSACVoltageSource): amplitude = bc.condition.source.amplitude return (self.analysis_spec.freqs, amplitude) + + def _get_normal_axis(self) -> int: + """Get the axis normal to the 2D simulation plane. + + Returns + ------- + int + Axis index (0=x, 1=y, 2=z) of the zero-size dimension. + + Raises + ------ + SetupError + If simulation is not 2D (exactly one zero-size dimension required). + """ + zero_dims = [i for i, s in enumerate(self.size) if s == 0] + if len(zero_dims) != 1: + raise SetupError( + "Automatic doping refinement is only supported for 2D simulations " + f"(exactly one zero-size dimension). Found {len(zero_dims)} zero dimensions." + ) + return zero_dims[0] + + def generate_doping_refinement_regions(self) -> list[GridRefinementRegion]: + """Generate mesh refinement regions based on doping gradients in semiconductor structures. + + This method analyzes all semiconductor structures in the simulation and generates + :class:`GridRefinementRegion` objects in areas where doping concentration varies + rapidly. The generated regions can be added to the grid specification. + + Returns + ------- + list[GridRefinementRegion] + List of refinement regions. Empty if no auto refinement is configured or + no semiconductor structures are present. + + Notes + ----- + This method requires: + - A 2D simulation (exactly one zero-size dimension) + - ``DistanceUnstructuredGrid`` with ``auto_doping_refinement`` specified + - At least one structure with ``SemiconductorMedium`` + + Example + ------- + >>> # Get auto-generated refinement regions + >>> regions = sim.generate_doping_refinement_regions() + >>> # Add to existing mesh refinements + >>> new_grid = sim.grid_spec.updated_copy( + ... mesh_refinements=sim.grid_spec.mesh_refinements + tuple(regions) + ... ) + """ + refinement_regions = [] + + # Check if auto refinement is configured + if not isinstance(self.grid_spec, DistanceUnstructuredGrid): + return refinement_regions + + spec = self.grid_spec.auto_doping_refinement + if spec is None: + return refinement_regions + + # Check for 2D simulation + try: + normal_axis = self._get_normal_axis() + except SetupError: + log.warning( + "Automatic doping refinement is only supported for 2D simulations. " + "Skipping auto refinement." + ) + return refinement_regions + + # Process each semiconductor structure + for structure in self.structures: + medium = structure.medium + + # Handle MultiPhysicsMedium + if isinstance(medium, MultiPhysicsMedium): + charge_medium = medium.charge + else: + charge_medium = medium if isinstance(medium, SemiconductorMedium) else None + + if not isinstance(charge_medium, SemiconductorMedium): + continue + + # Check if there's non-trivial doping + has_doping = (isinstance(charge_medium.N_d, tuple) and len(charge_medium.N_d) > 0) or ( + isinstance(charge_medium.N_a, tuple) and len(charge_medium.N_a) > 0 + ) + + if not has_doping: + continue + + # Get structure bounds, clipped to simulation domain + struct_bounds = structure.geometry.bounds + sim_bounds = self.bounds + + clipped_bounds = ( + tuple(max(struct_bounds[0][i], sim_bounds[0][i]) for i in range(3)), + tuple(min(struct_bounds[1][i], sim_bounds[1][i]) for i in range(3)), + ) + + # Generate refinement regions for this structure + struct_regions = charge_medium.compute_doping_refinement_regions( + bounds=clipped_bounds, + spec=spec, + normal_axis=normal_axis, + ) + + refinement_regions.extend(struct_regions) + + return refinement_regions + + def with_auto_doping_refinement(self) -> HeatChargeSimulation: + """Return a copy of this simulation with auto-generated doping refinement regions. + + This is a convenience method that generates refinement regions based on + doping gradients and returns an updated simulation with these regions + added to the grid specification. + + Returns + ------- + HeatChargeSimulation + Copy of this simulation with additional refinement regions. + + Example + ------- + >>> sim_refined = sim.with_auto_doping_refinement() + """ + if not isinstance(self.grid_spec, DistanceUnstructuredGrid): + return self + + new_regions = self.generate_doping_refinement_regions() + + if not new_regions: + return self + + combined_refinements = list(self.grid_spec.mesh_refinements) + new_regions + + # Merge overlapping regions by taking minimum dl_internal + merged_refinements = self._merge_refinement_regions(combined_refinements) + + new_grid_spec = self.grid_spec.updated_copy(mesh_refinements=tuple(merged_refinements)) + + return self.updated_copy(grid_spec=new_grid_spec) + + @staticmethod + def _merge_refinement_regions( + regions: list[GridRefinementRegion], + ) -> list[GridRefinementRegion]: + """Merge overlapping refinement regions by taking minimum dl_internal. + + Currently implements a simple approach: no merging, just return as-is. + Future enhancement could implement proper spatial merging. + """ + # For now, just return the regions as-is + # A more sophisticated approach would merge overlapping boxes + return regions From 8853fe6cd33c2a55a2c315ac9cdbf5db71a07e32 Mon Sep 17 00:00:00 2001 From: Momchil Minkov Date: Mon, 15 Dec 2025 13:54:05 +0100 Subject: [PATCH 2/3] fixes --- docs/notebooks | 2 +- tests/test_components/test_heat_charge.py | 219 ++++++++++++++++++++++ tidy3d/components/material/tcad/charge.py | 12 +- tidy3d/components/tcad/doping.py | 21 +++ 4 files changed, 247 insertions(+), 7 deletions(-) diff --git a/docs/notebooks b/docs/notebooks index c9a54bcb87..2571f1e39e 160000 --- a/docs/notebooks +++ b/docs/notebooks @@ -1 +1 @@ -Subproject commit c9a54bcb87a25de06bda0c41f22d84ffd8605cd3 +Subproject commit 2571f1e39e95a7a894beb655884ef869a564dc48 diff --git a/tests/test_components/test_heat_charge.py b/tests/test_components/test_heat_charge.py index 6d0ed00655..2c2ba1cd9d 100644 --- a/tests/test_components/test_heat_charge.py +++ b/tests/test_components/test_heat_charge.py @@ -2597,3 +2597,222 @@ def test_heat_only_simulation_with_semiconductor(): assert TCADAnalysisTypes.CONDUCTION not in simulation_types, ( "Conduction simulation should NOT be triggered when no electric BCs are present." ) + + +def test_auto_doping_refinement(): + """Test automatic doping-based mesh refinement generation.""" + # Create a semiconductor medium with Gaussian doping (pn junction) + Si_doped = td.MultiPhysicsMedium( + charge=td.SemiconductorMedium( + permittivity=11.7, + N_d=( + td.GaussianDoping( + center=(0.1, 0, 0), + size=(0.2, 0.5, td.inf), # Use inf for the zero dimension + ref_con=1e15, + concentration=1e18, + width=0.05, + source="xmin", + ), + ), + N_a=( + td.GaussianDoping( + center=(-0.1, 0, 0), + size=(0.2, 0.5, td.inf), # Use inf for the zero dimension + ref_con=1e15, + concentration=1e18, + width=0.05, + source="xmax", + ), + ), + N_c=td.ConstantEffectiveDOS(N=2.86e19), + N_v=td.ConstantEffectiveDOS(N=3.1e19), + E_g=td.ConstantEnergyBandGap(eg=1.11), + mobility_n=td.ConstantMobilityModel(mu=1500), + mobility_p=td.ConstantMobilityModel(mu=500), + ), + name="Si_doped", + ) + + # Create structure + semiconductor_structure = td.Structure( + geometry=td.Box(center=(0, 0, 0), size=(0.5, 0.5, td.inf)), + medium=Si_doped, + name="semiconductor", + ) + + # Create oxide medium (insulator) + oxide = td.MultiPhysicsMedium( + charge=td.ChargeInsulatorMedium(permittivity=3.9), + name="oxide", + ) + + # Create a grid spec with auto doping refinement enabled + grid_with_auto = td.DistanceUnstructuredGrid( + dl_interface=0.01, + dl_bulk=0.1, + distance_interface=0.02, + distance_bulk=0.2, + auto_doping_refinement=td.DopingGradientRefinementSpec( + resolution_factor=5.0, + min_dl=0.005, + max_dl=0.05, + num_samples=30, + transition_factor=10.0, + ), + ) + + # Create a 2D HeatChargeSimulation + sim = td.HeatChargeSimulation( + medium=oxide, + structures=[semiconductor_structure], + center=(0, 0, 0), + size=(1, 1, 0), # 2D simulation (z=0) + boundary_spec=[], + grid_spec=grid_with_auto, + sources=[], + monitors=[], + ) + + # Test generate_doping_refinement_regions + regions = sim.generate_doping_refinement_regions() + + # Should generate some refinement regions for the doping gradients + assert isinstance(regions, list), "Should return a list of regions" + + # Check that regions are GridRefinementRegion objects + for region in regions: + assert isinstance(region, td.GridRefinementRegion), ( + f"Expected GridRefinementRegion, got {type(region)}" + ) + # Check that dl_internal is within expected bounds + assert 0.005 <= region.dl_internal <= 0.05, ( + f"dl_internal {region.dl_internal} outside expected range [0.005, 0.05]" + ) + + # Test with_auto_doping_refinement + sim_refined = sim.with_auto_doping_refinement() + assert isinstance(sim_refined, td.HeatChargeSimulation), "Should return a HeatChargeSimulation" + + # The refined simulation should have the auto-generated regions added to mesh_refinements + if len(regions) > 0: + assert len(sim_refined.grid_spec.mesh_refinements) >= len(regions), ( + "Refined simulation should have at least as many mesh refinements as auto-generated regions" + ) + + +def test_auto_doping_refinement_disabled(): + """Test that auto doping refinement can be disabled.""" + Si_doped = td.MultiPhysicsMedium( + charge=td.SemiconductorMedium( + permittivity=11.7, + N_d=( + td.GaussianDoping( + center=(0.1, 0, 0), + size=(0.2, 0.5, td.inf), # Use inf for the zero dimension + ref_con=1e15, + concentration=1e18, + width=0.05, + source="xmin", + ), + ), + N_c=td.ConstantEffectiveDOS(N=2.86e19), + N_v=td.ConstantEffectiveDOS(N=3.1e19), + E_g=td.ConstantEnergyBandGap(eg=1.11), + mobility_n=td.ConstantMobilityModel(mu=1500), + mobility_p=td.ConstantMobilityModel(mu=500), + ), + name="Si_doped", + ) + + semiconductor_structure = td.Structure( + geometry=td.Box(center=(0, 0, 0), size=(0.5, 0.5, td.inf)), + medium=Si_doped, + name="semiconductor", + ) + + oxide = td.MultiPhysicsMedium( + charge=td.ChargeInsulatorMedium(permittivity=3.9), + name="oxide", + ) + + # Grid with auto_doping_refinement explicitly disabled + grid_no_auto = td.DistanceUnstructuredGrid( + dl_interface=0.01, + dl_bulk=0.1, + distance_interface=0.02, + distance_bulk=0.2, + auto_doping_refinement=None, + ) + + sim = td.HeatChargeSimulation( + medium=oxide, + structures=[semiconductor_structure], + center=(0, 0, 0), + size=(1, 1, 0), + boundary_spec=[], + grid_spec=grid_no_auto, + sources=[], + monitors=[], + ) + + # Should return empty list when auto_doping_refinement is None + regions = sim.generate_doping_refinement_regions() + assert regions == [], "Should return empty list when auto_doping_refinement is None" + + +def test_auto_doping_refinement_3d_warning(): + """Test that auto doping refinement warns for 3D simulations.""" + Si_doped = td.MultiPhysicsMedium( + charge=td.SemiconductorMedium( + permittivity=11.7, + N_d=( + td.ConstantDoping( + center=(0.1, 0, 0), + size=(0.2, 0.5, 0.2), + concentration=1e18, + ), + ), + N_c=td.ConstantEffectiveDOS(N=2.86e19), + N_v=td.ConstantEffectiveDOS(N=3.1e19), + E_g=td.ConstantEnergyBandGap(eg=1.11), + mobility_n=td.ConstantMobilityModel(mu=1500), + mobility_p=td.ConstantMobilityModel(mu=500), + ), + name="Si_doped", + ) + + semiconductor_structure = td.Structure( + geometry=td.Box(center=(0, 0, 0), size=(0.5, 0.5, 0.5)), + medium=Si_doped, + name="semiconductor", + ) + + oxide = td.MultiPhysicsMedium( + charge=td.ChargeInsulatorMedium(permittivity=3.9), + name="oxide", + ) + + grid_with_auto = td.DistanceUnstructuredGrid( + dl_interface=0.01, + dl_bulk=0.1, + distance_interface=0.02, + distance_bulk=0.2, + auto_doping_refinement=td.DopingGradientRefinementSpec(), + ) + + # 3D simulation + sim = td.HeatChargeSimulation( + medium=oxide, + structures=[semiconductor_structure], + center=(0, 0, 0), + size=(1, 1, 1), # 3D simulation + boundary_spec=[], + grid_spec=grid_with_auto, + sources=[], + monitors=[], + ) + + # Should return empty list for 3D simulation (with warning logged) + regions = sim.generate_doping_refinement_regions() + assert regions == [], "Should return empty list for 3D simulations" diff --git a/tidy3d/components/material/tcad/charge.py b/tidy3d/components/material/tcad/charge.py index fec0e0bf8f..45ae3b2653 100644 --- a/tidy3d/components/material/tcad/charge.py +++ b/tidy3d/components/material/tcad/charge.py @@ -464,13 +464,13 @@ def compute_doping_refinement_regions( grad_magnitude = grad_magnitude.squeeze() # Compute local mesh size: dl = N / (|grad_N| * k) - # Avoid division by zero + # Avoid division by zero by setting gradient to a small value where it's near zero eps = 1e-10 - local_dl = np.where( - grad_magnitude > eps, - net_doping / (grad_magnitude * spec.resolution_factor), - spec.max_dl, - ) + grad_safe = np.maximum(grad_magnitude, eps) + local_dl = net_doping / (grad_safe * spec.resolution_factor) + + # Where gradient was near zero, use max_dl + local_dl = np.where(grad_magnitude > eps, local_dl, spec.max_dl) # Clamp to [min_dl, max_dl] local_dl = np.clip(local_dl, spec.min_dl, spec.max_dl) diff --git a/tidy3d/components/tcad/doping.py b/tidy3d/components/tcad/doping.py index 439ba54f30..ee8754ca99 100644 --- a/tidy3d/components/tcad/doping.py +++ b/tidy3d/components/tcad/doping.py @@ -302,9 +302,20 @@ def gradient_magnitude(self, coords: dict, meshgrid: bool = True) -> np.ndarray: form involves products of Gaussians in each dimension. """ # Use finite differences on the doping profile + # Note: _get_contrib calls squeeze(), so we need to handle different dimensionalities doping = self._get_contrib(coords, meshgrid=meshgrid) if meshgrid: + # Ensure doping is 3D for consistent gradient computation + if doping.ndim == 2: + # Determine which axis was squeezed (has length 1 in coords) + if len(coords["x"]) == 1: + doping = doping[np.newaxis, :, :] + elif len(coords["y"]) == 1: + doping = doping[:, np.newaxis, :] + else: + doping = doping[:, :, np.newaxis] + # Compute gradient using finite differences dx = coords["x"][1] - coords["x"][0] if len(coords["x"]) > 1 else 1.0 dy = coords["y"][1] - coords["y"][0] if len(coords["y"]) > 1 else 1.0 @@ -399,6 +410,16 @@ def gradient_magnitude(self, coords: dict, meshgrid: bool = True) -> np.ndarray: doping = self._get_contrib(coords, meshgrid=meshgrid) if meshgrid: + # Ensure doping is 3D for consistent gradient computation + if doping.ndim == 2: + # Determine which axis was squeezed (has length 1 in coords) + if len(coords["x"]) == 1: + doping = doping[np.newaxis, :, :] + elif len(coords["y"]) == 1: + doping = doping[:, np.newaxis, :] + else: + doping = doping[:, :, np.newaxis] + dx = coords["x"][1] - coords["x"][0] if len(coords["x"]) > 1 else 1.0 dy = coords["y"][1] - coords["y"][0] if len(coords["y"]) > 1 else 1.0 dz = coords["z"][1] - coords["z"][0] if len(coords["z"]) > 1 else 1.0 From 490cad9331e5df83a9f273149ee36aab2d92faf5 Mon Sep 17 00:00:00 2001 From: Momchil Minkov Date: Mon, 15 Dec 2025 15:44:26 +0100 Subject: [PATCH 3/3] Add semiconductor interface grid refinement option --- docs/notebooks | 2 +- tests/test_components/test_heat_charge.py | 33 +++++++++++++++++++++++ tidy3d/components/tcad/grid.py | 10 +++++++ 3 files changed, 44 insertions(+), 1 deletion(-) diff --git a/docs/notebooks b/docs/notebooks index 2571f1e39e..c0a05cb924 160000 --- a/docs/notebooks +++ b/docs/notebooks @@ -1 +1 @@ -Subproject commit 2571f1e39e95a7a894beb655884ef869a564dc48 +Subproject commit c0a05cb92439038266a53c36e034ae18eb6e0adb diff --git a/tests/test_components/test_heat_charge.py b/tests/test_components/test_heat_charge.py index 2c2ba1cd9d..b8e7b309ed 100644 --- a/tests/test_components/test_heat_charge.py +++ b/tests/test_components/test_heat_charge.py @@ -2816,3 +2816,36 @@ def test_auto_doping_refinement_3d_warning(): # Should return empty list for 3D simulation (with warning logged) regions = sim.generate_doping_refinement_regions() assert regions == [], "Should return empty list for 3D simulations" + + +def test_dl_interface_semiconductor(): + """Test that dl_interface_semiconductor parameter works correctly.""" + + # Test default behavior (None defaults to dl_interface) + grid1 = td.DistanceUnstructuredGrid( + dl_interface=0.1, + dl_bulk=1.0, + distance_interface=0.3, + distance_bulk=2.0, + ) + assert grid1.dl_interface_semiconductor is None + + # Test explicit setting + grid2 = td.DistanceUnstructuredGrid( + dl_interface=0.1, + dl_interface_semiconductor=0.02, + dl_bulk=1.0, + distance_interface=0.3, + distance_bulk=2.0, + ) + assert grid2.dl_interface_semiconductor == 0.02 + + # Test that a smaller value can be set + grid3 = td.DistanceUnstructuredGrid( + dl_interface=0.1, + dl_interface_semiconductor=0.05, + dl_bulk=1.0, + distance_interface=0.3, + distance_bulk=2.0, + ) + assert grid3.dl_interface_semiconductor < grid3.dl_interface diff --git a/tidy3d/components/tcad/grid.py b/tidy3d/components/tcad/grid.py index 8a94d01e89..0171923e43 100644 --- a/tidy3d/components/tcad/grid.py +++ b/tidy3d/components/tcad/grid.py @@ -241,6 +241,16 @@ class DistanceUnstructuredGrid(UnstructuredGrid): units=MICROMETER, ) + dl_interface_semiconductor: pd.PositiveFloat = pd.Field( + None, + title="Semiconductor Interface Grid Size", + description="Grid size near interfaces between semiconductors and insulators (e.g., " + "Si/SiO2 interfaces). If ``None``, defaults to ``dl_interface``. Typically this should " + "be set to a smaller value than ``dl_interface`` to accurately resolve depletion " + "regions and charge accumulation at semiconductor-insulator boundaries.", + units=MICROMETER, + ) + dl_bulk: pd.PositiveFloat = pd.Field( ..., title="Bulk Grid Size",