""" Defines classes specifying meshing in 1D and a collective class for 3D """from__future__importannotationsfromabcimportABC,abstractmethodfromtypingimportTuple,List,Unionimportnumpyasnpimportpydantic.v1aspdfrom.gridimportCoords1D,Coords,Gridfrom.mesherimportGradedMesher,MesherTypefrom..baseimportTidy3dBaseModel,cached_propertyfrom..typesimportAxis,Symmetry,annotate_type,TYPE_TAG_STRfrom..sourceimportSourceTypefrom..structureimportStructure,StructureTypefrom..geometry.baseimportBoxfrom...logimportlogfrom...exceptionsimportSetupErrorfrom...constantsimportMICROMETER,C_0,fp_eps# Default Courant number reduction rate in Benkler's schemeDEFAULT_COURANT_REDUCTION_BENKLER=0.3classConformalMeshSpec(Tidy3dBaseModel,ABC):"""Base class defining conformal mesh specifications."""@cached_propertydefcourant_ratio(self)->float:"""The scaling ratio applied to Courant number so that the courant number in the simulation is ``sim.courant * courant_ratio``. """return1.0
[docs]classStaircasingConformalMeshSpec(ConformalMeshSpec):"""Simple staircasing scheme based on [Taflove, The Electrical Engineering Handbook 3.629-670 (2005): 15.]. """
[docs]classHeuristicConformalMeshSpec(ConformalMeshSpec):"""Slightly different from the staircasing scheme: the field component near PEC is considered to be outside PEC if it's substantially normal to the interface. """
[docs]classBenklerConformalMeshSpec(ConformalMeshSpec):"""Conformal mesh scheme based on [S. Benkler, IEEE Transactions on Antennas and Propagation 54.6, 1843 (2006)], which is similar to the approach described in [S. Dey, R. Mittra, IEEE Microwave and Guided Wave Letters 7.9, 273 (1997)]. """timestep_reduction:float=pd.Field(DEFAULT_COURANT_REDUCTION_BENKLER,title="Time Step Size Reduction Rate",description="Reduction factor between 0 and 1 such that the simulation's time step size ""will be ``1 - timestep_reduction`` times its default value. ""Accuracy can be improved with a smaller time step size; but simulation time increased as well.",lt=1,ge=0,)@cached_propertydefcourant_ratio(self)->float:"""The scaling ratio applied to Courant number so that the courant number in the simulation is ``sim.courant * courant_ratio``. """return1-self.timestep_reduction
[docs]classGridSpec1d(Tidy3dBaseModel,ABC):"""Abstract base class, defines 1D grid generation specifications."""
[docs]defmake_coords(self,axis:Axis,structures:List[StructureType],symmetry:Tuple[Symmetry,Symmetry,Symmetry],periodic:bool,wavelength:pd.PositiveFloat,num_pml_layers:Tuple[pd.NonNegativeInt,pd.NonNegativeInt],)->Coords1D:"""Generate 1D coords to be used as grid boundaries, based on simulation parameters. Symmetry, and PML layers will be treated here. Parameters ---------- axis : Axis Axis of this direction. structures : List[StructureType] List of structures present in simulation, the first one being the simulation domain. symmetry : Tuple[Symmetry, Symmetry, Symmetry] Reflection symmetry across a plane bisecting the simulation domain normal to each of the three axes. wavelength : float Free-space wavelength. num_pml_layers : Tuple[int, int] number of layers in the absorber + and - direction along one dimension. Returns ------- :class:`.Coords1D`: 1D coords to be used as grid boundaries. """# Determine if one should apply periodic boundary condition.# This should only affect auto nonuniform mesh generation for now.is_periodic=periodicandsymmetry[axis]==0# generate boundariesbound_coords=self._make_coords_initial(axis=axis,structures=structures,wavelength=wavelength,symmetry=symmetry,is_periodic=is_periodic,)# incorporate symmetriesifsymmetry[axis]!=0:# Offset to center if symmetry presentcenter=structures[0].geometry.center[axis]center_ind=np.argmin(np.abs(center-bound_coords))bound_coords+=center-bound_coords[center_ind]bound_coords=bound_coords[bound_coords>=center]bound_coords=np.append(2*center-bound_coords[:0:-1],bound_coords)# Add PML layers in using dl on edgesbound_coords=self._add_pml_to_bounds(num_pml_layers,bound_coords)returnbound_coords
@abstractmethoddef_make_coords_initial(self,axis:Axis,structures:List[StructureType],**kwargs,)->Coords1D:"""Generate 1D coords to be used as grid boundaries, based on simulation parameters. Symmetry, PML etc. are not considered in this method. For auto nonuniform generation, it will take some more arguments. Parameters ---------- structures : List[StructureType] List of structures present in simulation, the first one being the simulation domain. **kwargs Other arguments Returns ------- :class:`.Coords1D`: 1D coords to be used as grid boundaries. """@staticmethoddef_add_pml_to_bounds(num_layers:Tuple[int,int],bounds:Coords1D)->Coords1D:"""Append absorber layers to the beginning and end of the simulation bounds along one dimension. Parameters ---------- num_layers : Tuple[int, int] number of layers in the absorber + and - direction along one dimension. bound_coords : np.ndarray coordinates specifying boundaries between cells along one dimension. Returns ------- np.ndarray New bound coordinates along dimension taking abosrber into account. """ifbounds.size<2:returnboundsfirst_step=bounds[1]-bounds[0]last_step=bounds[-1]-bounds[-2]add_left=bounds[0]-first_step*np.arange(num_layers[0],0,-1)add_right=bounds[-1]+last_step*np.arange(1,num_layers[1]+1)returnnp.concatenate((add_left,bounds,add_right))
[docs]classUniformGrid(GridSpec1d):"""Uniform 1D grid. The most standard way to define a simulation is to use a constant grid size in each of the three directions. Example ------- >>> grid_1d = UniformGrid(dl=0.1) See Also -------- :class:`AutoGrid` Specification for non-uniform grid along a given dimension. **Notebooks:** * `Photonic crystal waveguide polarization filter <../../../notebooks/PhotonicCrystalWaveguidePolarizationFilter.html>`_ * `Using automatic nonuniform meshing <../../notebooks/AutoGrid.html>`_ """dl:pd.PositiveFloat=pd.Field(...,title="Grid Size",description="Grid size for uniform grid generation.",units=MICROMETER,)def_make_coords_initial(self,axis:Axis,structures:List[StructureType],**kwargs,)->Coords1D:"""Uniform 1D coords to be used as grid boundaries. Parameters ---------- axis : Axis Axis of this direction. structures : List[StructureType] List of structures present in simulation, the first one being the simulation domain. **kwargs: Other arguments all go here. Returns ------- :class:`.Coords1D`: 1D coords to be used as grid boundaries. """center,size=structures[0].geometry.center[axis],structures[0].geometry.size[axis]# Take a number of steps commensurate with the size; make dl a bit smaller if needednum_cells=int(np.ceil(size/self.dl))# Make sure there's at least one cellnum_cells=max(num_cells,1)# Adjust step size to fit simulation size exactlydl_snapped=size/num_cellsifsize>0elseself.dlreturncenter-size/2+np.arange(num_cells+1)*dl_snapped
[docs]classCustomGrid(GridSpec1d):"""Custom 1D grid supplied as a list of grid cell sizes centered on the simulation center. Example ------- >>> grid_1d = CustomGrid(dl=[0.2, 0.2, 0.1, 0.1, 0.1, 0.2, 0.2]) """dl:Tuple[pd.PositiveFloat,...]=pd.Field(...,title="Customized grid sizes.",description="An array of custom nonuniform grid sizes. The resulting grid is centered on ""the simulation center such that it spans the region ""``(center - sum(dl)/2, center + sum(dl)/2)``, unless a ``custom_offset`` is given. ""Note: if supplied sizes do not cover the simulation size, the first and last sizes ""are repeated to cover the simulation domain.",units=MICROMETER,)custom_offset:float=pd.Field(None,title="Customized grid offset.",description="The starting coordinate of the grid which defines the simulation center. ""If ``None``, the simulation center is set such that it spans the region ""``(center - sum(dl)/2, center + sum(dl)/2)``.",units=MICROMETER,)def_make_coords_initial(self,axis:Axis,structures:List[StructureType],**kwargs,)->Coords1D:"""Customized 1D coords to be used as grid boundaries. Parameters ---------- axis : Axis Axis of this direction. structures : List[StructureType] List of structures present in simulation, the first one being the simulation domain. *kwargs Other arguments all go here. Returns ------- :class:`.Coords1D`: 1D coords to be used as grid boundaries. """center,size=structures[0].geometry.center[axis],structures[0].geometry.size[axis]# get bounding coordinatesdl=np.array(self.dl)bound_coords=np.append(0.0,np.cumsum(dl))# place the middle of the bounds at the center of the simulation along dimension,# or use the `custom_offset` if providedifself.custom_offsetisNone:bound_coords+=center-bound_coords[-1]/2else:bound_coords+=self.custom_offset# chop off any coords outside of simulation bounds, beyond some buffer region# to take numerical effects into accountbuffer=fp_eps*sizebound_min=center-size/2-bufferbound_max=center+size/2+bufferifbound_max<bound_coords[0]orbound_min>bound_coords[-1]:axis_name="xyz"[axis]raiseSetupError(f"Simulation domain does not overlap with the provided custom grid in '{axis_name}' direction.")ifsize==0:# in case of zero-size dimension return the boundaries between which simulation fallsind=np.searchsorted(bound_coords,center,side="right")# in case when the center coincides with the right most boundaryifind>=len(bound_coords):ind=len(bound_coords)-1returnbound_coords[ind-1:ind+1]else:bound_coords=bound_coords[bound_coords<=bound_max]bound_coords=bound_coords[bound_coords>=bound_min]# if not extending to simulation bounds, repeat beginning and enddl_min=dl[0]dl_max=dl[-1]whilebound_coords[0]-dl_min>=bound_min:bound_coords=np.insert(bound_coords,0,bound_coords[0]-dl_min)whilebound_coords[-1]+dl_max<=bound_max:bound_coords=np.append(bound_coords,bound_coords[-1]+dl_max)# in case a `custom_offset` is provided, it's possible the bounds were numerically within# the simulation bounds but were still chopped off, which is fixed hereifself.custom_offsetisnotNone:ifnp.isclose(bound_coords[0]-dl_min,bound_min):bound_coords=np.insert(bound_coords,0,bound_coords[0]-dl_min)ifnp.isclose(bound_coords[-1]+dl_max,bound_max):bound_coords=np.append(bound_coords,bound_coords[-1]+dl_max)returnbound_coords
[docs]classAutoGrid(GridSpec1d):"""Specification for non-uniform grid along a given dimension. Example ------- >>> grid_1d = AutoGrid(min_steps_per_wvl=16, max_scale=1.4) See Also -------- :class:`UniformGrid` Uniform 1D grid. :class:`GridSpec` Collective grid specification for all three dimensions. **Notebooks:** * `Using automatic nonuniform meshing <../../notebooks/AutoGrid.html>`_ **Lectures:** * `Time step size and CFL condition in FDTD <https://www.flexcompute.com/fdtd101/Lecture-7-Time-step-size-and-CFL-condition-in-FDTD/>`_ * `Numerical dispersion in FDTD <https://www.flexcompute.com/fdtd101/Lecture-8-Numerical-dispersion-in-FDTD/>`_ """min_steps_per_wvl:float=pd.Field(10.0,title="Minimal number of steps per wavelength",description="Minimal number of steps per wavelength in each medium.",ge=6.0,)max_scale:float=pd.Field(1.4,title="Maximum Grid Size Scaling",description="Sets the maximum ratio between any two consecutive grid steps.",ge=1.2,lt=2.0,)dl_min:pd.NonNegativeFloat=pd.Field(0,title="Lower bound of grid size",description="Lower bound of the grid size along this dimension regardless of ""structures present in the simulation, including override structures ""with ``enforced=True``. It is a soft bound, meaning that the actual minimal ""grid size might be slightly smaller.",units=MICROMETER,)mesher:MesherType=pd.Field(GradedMesher(),title="Grid Construction Tool",description="The type of mesher to use to generate the grid automatically.",)def_make_coords_initial(self,axis:Axis,structures:List[StructureType],wavelength:float,symmetry:Symmetry,is_periodic:bool,)->Coords1D:"""Customized 1D coords to be used as grid boundaries. Parameters ---------- axis : Axis Axis of this direction. structures : List[StructureType] List of structures present in simulation. wavelength : float Free-space wavelength. symmetry : Tuple[Symmetry, Symmetry, Symmetry] Reflection symmetry across a plane bisecting the simulation domain normal to each of the three axes. is_periodic : bool Apply periodic boundary condition or not. Returns ------- :class:`.Coords1D`: 1D coords to be used as grid boundaries. """sim_cent=list(structures[0].geometry.center)sim_size=list(structures[0].geometry.size)fordim,syminenumerate(symmetry):ifsym!=0:sim_cent[dim]+=sim_size[dim]/4sim_size[dim]/=2symmetry_domain=Box(center=sim_cent,size=sim_size)# New list of structures with symmetry appliedstruct_list=[Structure(geometry=symmetry_domain,medium=structures[0].medium)]forstructureinstructures[1:]:ifsymmetry_domain.intersects(structure.geometry):struct_list.append(structure)# parse structuresinterval_coords,max_dl_list=self.mesher.parse_structures(axis,struct_list,wavelength,self.min_steps_per_wvl,self.dl_min,)# Put just a single pixel if 2D-like simulationifinterval_coords.size==1:dl=wavelength/self.min_steps_per_wvlreturnnp.array([sim_cent[axis]-dl/2,sim_cent[axis]+dl/2])# generate mesh stepsinterval_coords=np.array(interval_coords).flatten()max_dl_list=np.array(max_dl_list).flatten()len_interval_list=interval_coords[1:]-interval_coords[:-1]dl_list=self.mesher.make_grid_multiple_intervals(max_dl_list,len_interval_list,self.max_scale,is_periodic)# generate boundariesbound_coords=np.append(0.0,np.cumsum(np.concatenate(dl_list)))bound_coords+=interval_coords[0]# fix simulation domain boundaries which may be slightly offdomain_bounds=[bound[axis]forboundinsymmetry_domain.bounds]ifnotnp.all(np.isclose(bound_coords[[0,-1]],domain_bounds)):raiseSetupError(f"AutoGrid coordinates along axis {axis} do not match the simulation ""domain, indicating an unexpected error in the meshing. Please create a github ""issue so that the problem can be investigated. In the meantime, switch to "f"uniform or custom grid along {axis}.")bound_coords[[0,-1]]=domain_boundsreturnnp.array(bound_coords)
GridType=Union[UniformGrid,CustomGrid,AutoGrid]
[docs]classGridSpec(Tidy3dBaseModel):"""Collective grid specification for all three dimensions. Example ------- >>> uniform = UniformGrid(dl=0.1) >>> custom = CustomGrid(dl=[0.2, 0.2, 0.1, 0.1, 0.1, 0.2, 0.2]) >>> auto = AutoGrid(min_steps_per_wvl=12) >>> grid_spec = GridSpec(grid_x=uniform, grid_y=custom, grid_z=auto, wavelength=1.5) See Also -------- :class:`UniformGrid` Uniform 1D grid. :class:`AutoGrid` Specification for non-uniform grid along a given dimension. **Notebooks:** * `Using automatic nonuniform meshing <../../notebooks/AutoGrid.html>`_ **Lectures:** * `Time step size and CFL condition in FDTD <https://www.flexcompute.com/fdtd101/Lecture-7-Time-step-size-and-CFL-condition-in-FDTD/>`_ * `Numerical dispersion in FDTD <https://www.flexcompute.com/fdtd101/Lecture-8-Numerical-dispersion-in-FDTD/>`_ """grid_x:GridType=pd.Field(AutoGrid(),title="Grid specification along x-axis",description="Grid specification along x-axis",discriminator=TYPE_TAG_STR,)grid_y:GridType=pd.Field(AutoGrid(),title="Grid specification along y-axis",description="Grid specification along y-axis",discriminator=TYPE_TAG_STR,)grid_z:GridType=pd.Field(AutoGrid(),title="Grid specification along z-axis",description="Grid specification along z-axis",discriminator=TYPE_TAG_STR,)wavelength:float=pd.Field(None,title="Free-space wavelength",description="Free-space wavelength for automatic nonuniform grid. It can be 'None' ""if there is at least one source in the simulation, in which case it is defined by ""the source central frequency. ""Note: it only takes effect when at least one of the three dimensions ""uses :class:`.AutoGrid`.",units=MICROMETER,)override_structures:Tuple[annotate_type(StructureType),...]=pd.Field((),title="Grid specification override structures",description="A set of structures that is added on top of the simulation structures in ""the process of generating the grid. This can be used to refine the grid or make it ""coarser depending than the expected need for higher/lower resolution regions. ""Note: it only takes effect when at least one of the three dimensions ""uses :class:`.AutoGrid`.",)@propertydefauto_grid_used(self)->bool:"""True if any of the three dimensions uses :class:`.AutoGrid`."""grid_list=[self.grid_x,self.grid_y,self.grid_z]returnnp.any([isinstance(mesh,AutoGrid)formeshingrid_list])@propertydefcustom_grid_used(self)->bool:"""True if any of the three dimensions uses :class:`.CustomGrid`."""grid_list=[self.grid_x,self.grid_y,self.grid_z]returnnp.any([isinstance(mesh,CustomGrid)formeshingrid_list])
[docs]@staticmethoddefwavelength_from_sources(sources:List[SourceType])->pd.PositiveFloat:"""Define a wavelength based on supplied sources. Called if auto mesh is used and ``self.wavelength is None``."""# no sourcesiflen(sources)==0:raiseSetupError("Automatic grid generation requires the input of 'wavelength' or sources.")# Use central frequency of sources, if any.freqs=np.array([source.source_time.freq0forsourceinsources])# multiple sources of different central frequenciesifnotnp.all(np.isclose(freqs,freqs[0])):raiseSetupError("Sources of different central frequencies are supplied. ""Please supply a 'wavelength' value for 'grid_spec'.")returnC_0/freqs[0]
@propertydefoverride_structures_used(self)->List[bool,bool,bool]:"""Along each axis, ``True`` if any override structure is used. However, it is still ``False`` if only :class:`.MeshOverrideStructure` is supplied, and their ``dl[axis]`` all take the ``None`` value. """# empty override_structure listiflen(self.override_structures)==0:return[False]*3override_used=[False]*3forstructureinself.override_structures:# override used in all axes if any `Structure` is presentifisinstance(structure,Structure):return[True]*3fordl_axis,dlinenumerate(structure.dl):if(notoverride_used[dl_axis])and(dlisnotNone):override_used[dl_axis]=Truereturnoverride_used
[docs]defmake_grid(self,structures:List[Structure],symmetry:Tuple[Symmetry,Symmetry,Symmetry],periodic:Tuple[bool,bool,bool],sources:List[SourceType],num_pml_layers:List[Tuple[pd.NonNegativeInt,pd.NonNegativeInt]],)->Grid:"""Make the entire simulation grid based on some simulation parameters. Parameters ---------- structures : List[Structure] List of structures present in the simulation. The first structure must be the simulation geometry with the simulation background medium. symmetry : Tuple[Symmetry, Symmetry, Symmetry] Reflection symmetry across a plane bisecting the simulation domain normal to each of the three axes. sources : List[SourceType] List of sources. num_pml_layers : List[Tuple[float, float]] List containing the number of absorber layers in - and + boundaries. Returns ------- Grid: Entire simulation grid. """# Set up wavelength for automatic mesh generation if needed.wavelength=self.wavelengthifwavelengthisNoneandself.auto_grid_used:wavelength=self.wavelength_from_sources(sources)log.info(f"Auto meshing using wavelength {wavelength:1.4f} defined from sources.")# Warn user if ``GridType`` along some axis is not ``AutoGrid`` and# ``override_structures`` is not empty. The override structures# are not effective along those axes.foraxis_ind,override_used_axis,grid_axisinzip(["x","y","z"],self.override_structures_used,[self.grid_x,self.grid_y,self.grid_z]):ifoverride_used_axisandnotisinstance(grid_axis,AutoGrid):log.warning(f"Override structures take no effect along {axis_ind}-axis. ""If intending to apply override structures to this axis, ""use 'AutoGrid'.",capture=False,)grids_1d=[self.grid_x,self.grid_y,self.grid_z]coords_dict={}foridim,(dim,grid_1d)inenumerate(zip("xyz",grids_1d)):coords_dict[dim]=grid_1d.make_coords(axis=idim,structures=list(structures)+list(self.override_structures),symmetry=symmetry,periodic=periodic[idim],wavelength=wavelength,num_pml_layers=num_pml_layers[idim],)coords=Coords(**coords_dict)returnGrid(boundaries=coords)
[docs]@classmethoddefauto(cls,wavelength:pd.PositiveFloat=None,min_steps_per_wvl:pd.PositiveFloat=10.0,max_scale:pd.PositiveFloat=1.4,override_structures:List[StructureType]=(),dl_min:pd.NonNegativeFloat=0.0,mesher:MesherType=GradedMesher(),)->GridSpec:"""Use the same :class:`AutoGrid` along each of the three directions. Parameters ---------- wavelength : pd.PositiveFloat, optional Free-space wavelength for automatic nonuniform grid. It can be 'None' if there is at least one source in the simulation, in which case it is defined by the source central frequency. min_steps_per_wvl : pd.PositiveFloat, optional Minimal number of steps per wavelength in each medium. max_scale : pd.PositiveFloat, optional Sets the maximum ratio between any two consecutive grid steps. override_structures : List[StructureType] A list of structures that is added on top of the simulation structures in the process of generating the grid. This can be used to refine the grid or make it coarser depending than the expected need for higher/lower resolution regions. dl_min: pd.NonNegativeFloat Lower bound of grid size. mesher : MesherType = GradedMesher() The type of mesher to use to generate the grid automatically. Returns ------- GridSpec :class:`GridSpec` with the same automatic nonuniform grid settings in each direction. """grid_1d=AutoGrid(min_steps_per_wvl=min_steps_per_wvl,max_scale=max_scale,dl_min=dl_min,mesher=mesher,)returncls(wavelength=wavelength,grid_x=grid_1d,grid_y=grid_1d,grid_z=grid_1d,override_structures=override_structures,)
[docs]@classmethoddefuniform(cls,dl:float)->GridSpec:"""Use the same :class:`UniformGrid` along each of the three directions. Parameters ---------- dl : float Grid size for uniform grid generation. Returns ------- GridSpec :class:`GridSpec` with the same uniform grid size in each direction. """grid_1d=UniformGrid(dl=dl)returncls(grid_x=grid_1d,grid_y=grid_1d,grid_z=grid_1d)