Source code for fdtd.sources

""" Sources are objects that inject the fields into the grid.

Available sources:

- PointSource
- LineSource

"""
## Imports

# other
from math import pi, sin

# typing
from .typing_ import Tuple, Number, ListOrSlice, List
from numpy import ndarray

# relatvie
from .grid import Grid
from .backend import backend as bd
from .waveforms import *
from .detectors import CurrentDetector

## PointSource class
[docs]class PointSource: """A source placed at a single point (grid cell) in the grid"""
[docs] def __init__( self, period: Number = 15, amplitude: float = 1.0, phase_shift: float = 0.0, name: str = None, pulse: bool = False, cycle: int = 5, hanning_dt: float = 10.0, ): """Create a LineSource with a gaussian profile Args: period: The period of the source. The period can be specified as integer [timesteps] or as float [seconds] amplitude: The electric field amplitude in simulation units phase_shift: The phase offset of the source. name: name of the source. pulse: Set True to use a Hanning window pulse instead of continuous wavefunction. cycle: cycles for Hanning window pulse. hanning_dt: timestep used for Hanning window pulse width (optional). """ self.grid = None self.period = period self.amplitude = amplitude self.phase_shift = phase_shift self.name = name self.pulse = pulse self.cycle = cycle self.frequency = 1.0 / period self.hanning_dt = hanning_dt if hanning_dt is not None else 0.5 / self.frequency
def _register_grid(self, grid: Grid, x: Number, y: Number, z: Number): """Register a grid for the source. Args: grid: the grid to place the source into. x: The x-location of the source in the grid y: The y-location of the source in the grid z: The z-location of the source in the grid Note: As its name suggests, this source is a POINT source. Hence it should be placed at a single coordinate tuple int the grid. """ self.grid = grid self.grid.sources.append(self) if self.name is not None: if not hasattr(grid, self.name): setattr(grid, self.name, self) else: raise ValueError( f"The grid already has an attribute with name {self.name}" ) try: (x,), (y,), (z,) = x, y, z except (TypeError, ValueError): raise ValueError("a point source should be placed on a single grid cell.") self.x, self.y, self.z = grid._handle_tuple((x, y, z)) self.period = grid._handle_time(self.period)
[docs] def update_E(self): """Add the source to the electric field""" q = self.grid.time_steps_passed # if pulse if self.pulse: t1 = int(2 * pi / (self.frequency * self.hanning_dt / self.cycle)) if q < t1: src = self.amplitude * hanning( self.frequency, q * self.hanning_dt, self.cycle ) else: # src = - self.grid.E[self.x, self.y, self.z, 2] src = 0 # if not pulse else: src = self.amplitude * sin(2 * pi * q / self.period + self.phase_shift) self.grid.E[self.x, self.y, self.z, 2] += src
[docs] def update_H(self): """Add the source to the magnetic field"""
def __repr__(self): return ( f"{self.__class__.__name__}(period={self.period}, " f"amplitude={self.amplitude}, phase_shift={self.phase_shift}, " f"name={repr(self.name)})" ) def __str__(self): s = " " + repr(self) + "\n" x = f"{self.x}" y = f"{self.y}" z = f"{self.z}" s += f" @ x={x}, y={y}, z={z}\n" return s
## LineSource class
[docs]class LineSource: """A source along a line in the FDTD grid"""
[docs] def __init__( self, period: Number = 15, amplitude: float = 1.0, phase_shift: float = 0.0, name: str = None, pulse: bool = False, cycle: int = 5, hanning_dt: float = 10.0, ): """Create a LineSource with a gaussian profile Args: period: The period of the source. The period can be specified as integer [timesteps] or as float [seconds] amplitude: The amplitude of the source in simulation units phase_shift: The phase offset of the source. pulse: Set True to use a Hanning window pulse instead of continuous wavefunction. cycle: cycles for Hanning window pulse. hanning_dt: timestep used for Hanning window pulse width (optional). """ self.grid = None self.period = period self.amplitude = amplitude self.phase_shift = phase_shift self.name = name self.pulse = pulse self.cycle = cycle self.frequency = 1.0 / period self.hanning_dt = hanning_dt if hanning_dt is not None else 0.5 / self.frequency
def _register_grid( self, grid: Grid, x: ListOrSlice, y: ListOrSlice, z: ListOrSlice ): """Register a grid for the source. Args: grid: the grid to place the source into. x: The x-location of the source in the grid y: The y-location of the source in the grid z: The z-location of the source in the grid Note: As its name suggests, this source is a LINE source. Hence the source spans the diagonal of the cube defined by the slices in the grid. """ self.grid = grid self.grid.sources.append(self) if self.name is not None: if not hasattr(grid, self.name): setattr(grid, self.name, self) else: raise ValueError( f"The grid already has an attribute with name {self.name}" ) self.x, self.y, self.z = self._handle_slices(x, y, z) self.period = grid._handle_time(self.period) L = len(self.x) vect = bd.array( (bd.array(self.x) - self.x[L // 2]) ** 2 + (bd.array(self.y) - self.y[L // 2]) ** 2 + (bd.array(self.z) - self.z[L // 2]) ** 2, bd.float, ) self.profile = bd.exp(-(vect ** 2) / (2 * (0.5 * vect.max()) ** 2)) self.profile /= self.profile.sum() self.profile *= self.amplitude def _handle_slices( self, x: ListOrSlice, y: ListOrSlice, z: ListOrSlice ) -> Tuple[List, List, List]: """Convert slices in the grid to lists This is necessary to make the source span the diagonal of the volume defined by the slices. Args: x: The x-location of the volume in the grid y: The y-location of the volume in the grid z: The z-location of the volume in the grid Returns: x, y, z: the x, y and z coordinates of the source as lists """ # if list-indices were chosen: if isinstance(x, list) and isinstance(y, list) and isinstance(z, list): if len(x) != len(y) or len(y) != len(z) or len(z) != len(x): raise IndexError( "sources require grid to be indexed with slices or equal length list-indices" ) x = [self.grid._handle_distance(_x) for _x in x] y = [self.grid._handle_distance(_y) for _y in y] z = [self.grid._handle_distance(_z) for _z in z] return x, y, z # if a combination of list-indices and slices were chosen, # convert the list-indices to slices. # TODO: maybe issue a warning here? if isinstance(x, list): x = slice( self.grid._handle_distance(x[0]), self.grid._handle_distance(x[-1]), None, ) if isinstance(y, list): y = slice( self.grid._handle_distance(y[0]), self.grid._handle_distance(y[-1]), None, ) if isinstance(z, list): z = slice( self.grid._handle_distance(z[0]), self.grid._handle_distance(z[-1]), None, ) # if we get here, we can assume slices: x0 = self.grid._handle_distance(x.start if x.start is not None else 0) y0 = self.grid._handle_distance(y.start if y.start is not None else 0) z0 = self.grid._handle_distance(z.start if z.start is not None else 0) x1 = self.grid._handle_distance(x.stop if x.stop is not None else self.grid.Nx) y1 = self.grid._handle_distance(y.stop if y.stop is not None else self.grid.Ny) z1 = self.grid._handle_distance(z.stop if z.stop is not None else self.grid.Nz) # we can now convert these coordinates into index lists m = max(abs(x1 - x0), abs(y1 - y0), abs(z1 - z0)) if m < 2: raise ValueError("a LineSource should consist of at least two gridpoints") x = [v.item() for v in bd.array(bd.linspace(x0, x1, m, endpoint=False), bd.int)] y = [v.item() for v in bd.array(bd.linspace(y0, y1, m, endpoint=False), bd.int)] z = [v.item() for v in bd.array(bd.linspace(z0, z1, m, endpoint=False), bd.int)] return x, y, z
[docs] def update_E(self): """Add the source to the electric field""" q = self.grid.time_steps_passed # if pulse if self.pulse: t1 = int(2 * pi / (self.frequency * self.hanning_dt / self.cycle)) if q < t1: vect = self.profile * hanning( self.frequency, q * self.hanning_dt, self.cycle ) else: # src = - self.grid.E[self.x, self.y, self.z, 2] vect = self.profile * 0 # if not pulse else: vect = self.profile * sin(2 * pi * q / self.period + self.phase_shift) # do not use list indexing here, as this is much slower especially for torch backend # DISABLED: self.grid.E[self.x, self.y, self.z, 2] = vect for x, y, z, value in zip(self.x, self.y, self.z, vect): self.grid.E[x, y, z, 2] += value
[docs] def update_H(self): """Add the source to the magnetic field"""
def __repr__(self): return ( f"{self.__class__.__name__}(period={self.period}, " f"amplitude={self.amplitude}, phase_shift={self.phase_shift}, " f"name={repr(self.name)})" ) def __str__(self): s = " " + repr(self) + "\n" x = f"[{self.x[0]}, ... , {self.x[-1]}]" y = f"[{self.y[0]}, ... , {self.y[-1]}]" z = f"[{self.z[0]}, ... , {self.z[-1]}]" s += f" @ x={x}, y={y}, z={z}\n" return s
## PlaneSource class
[docs]class PlaneSource: """A source along a plane in the FDTD grid"""
[docs] def __init__( self, period: Number = 15, amplitude: float = 1.0, phase_shift: float = 0.0, name: str = None, polarization: str = 'z', ): """Create a PlaneSource. Args: period: The period of the source. The period can be specified as integer [timesteps] or as float [seconds] amplitude: The amplitude of the source in simulation units phase_shift: The phase offset of the source. polarization: Axis of E-field polarization ('x','y',or 'z') """ self.grid = None self.period = period self.amplitude = amplitude self.phase_shift = phase_shift self.name = name self.polarization = polarization
def _register_grid( self, grid: Grid, x: ListOrSlice, y: ListOrSlice, z: ListOrSlice ): """Register a grid for the source. Args: grid: the grid to place the source into. x: The x-location of the source in the grid y: The y-location of the source in the grid z: The z-location of the source in the grid Note: As its name suggests, this source is a LINE source. Hence the source spans the diagonal of the cube defined by the slices in the grid. """ self.grid = grid self.grid.sources.append(self) if self.name is not None: if not hasattr(grid, self.name): setattr(grid, self.name, self) else: raise ValueError( f"The grid already has an attribute with name {self.name}" ) self.x, self.y, self.z = self._handle_slices(x, y, z) self.period = grid._handle_time(self.period) x = bd.arange(self.x.start, self.x.stop, 1) - (self.x.start + self.x.stop) // 2 y = bd.arange(self.y.start, self.y.stop, 1) - (self.y.start + self.y.stop) // 2 z = bd.arange(self.z.start, self.z.stop, 1) - (self.z.start + self.z.stop) // 2 xvec, yvec, zvec = bd.broadcast_arrays( x[:, None, None], y[None, :, None], z[None, None, :] ) _xvec = bd.array(xvec, float) _yvec = bd.array(yvec, float) _zvec = bd.array(zvec, float) profile = bd.ones(_xvec.shape) self.profile = self.amplitude * profile def _handle_slices( self, x: ListOrSlice, y: ListOrSlice, z: ListOrSlice ) -> Tuple[List, List, List]: """Validate slices and calculate center of plane Args: x: The x-location of the volume in the grid y: The y-location of the volume in the grid z: The z-location of the volume in the grid Returns: x, y, z: the x, y and z coordinates of the source as slices """ # ensure all slices if not isinstance(x, slice): if isinstance(x, list): (x,) = x x = slice( self.grid._handle_distance(x), self.grid._handle_distance(x) + 1, None ) if not isinstance(y, slice): if isinstance(y, list): (y,) = y y = slice( self.grid._handle_distance(y), self.grid._handle_distance(y) + 1, None ) if not isinstance(z, slice): if isinstance(z, list): (z,) = z z = slice( self.grid._handle_distance(z), self.grid._handle_distance(z) + 1, None ) # if we get here, we can assume slices: x0 = self.grid._handle_distance(x.start if x.start is not None else 0) y0 = self.grid._handle_distance(y.start if y.start is not None else 0) z0 = self.grid._handle_distance(z.start if z.start is not None else 0) x1 = self.grid._handle_distance(x.stop if x.stop is not None else self.grid.Nx) y1 = self.grid._handle_distance(y.stop if y.stop is not None else self.grid.Ny) z1 = self.grid._handle_distance(z.stop if z.stop is not None else self.grid.Nz) # make sure all slices have a start, stop and no step: x = ( slice(x0, x1) if x0 < x1 else (slice(x1, x0) if x0 > x1 else slice(x0, x0 + 1)) ) y = ( slice(y0, y1) if y0 < y1 else (slice(y1, y0) if y0 > y1 else slice(y0, y0 + 1)) ) z = ( slice(z0, z1) if z0 < z1 else (slice(z1, z0) if z0 > z1 else slice(z0, z0 + 1)) ) if [x.stop - x.start, y.stop - y.start, z.stop - z.start].count(0) > 0: raise ValueError( "Given location for PlaneSource results in slices of length 0!" ) if [x.stop - x.start, y.stop - y.start, z.stop - z.start].count(1) == 0: raise ValueError("Given location for PlaneSource is not a 2D plane!") if [x.stop - x.start, y.stop - y.start, z.stop - z.start].count(1) > 1: raise ValueError( "Given location for PlaneSource should have no more than one dimension in which it's flat.\n" "Use a LineSource for lower dimensional sources." ) self._Epol = 'xyz'.index(self.polarization) if (x.stop - x.start == 1 and self.polarization == 'x') or \ (y.stop - y.start == 1 and self.polarization == 'y') or \ (z.stop - z.start == 1 and self.polarization == 'z'): raise ValueError( "PlaneSource cannot be polarized perpendicular to the orientation of the plane." ) _Hpols = [(z,1,2), (z,0,2), (y,0,1)][self._Epol] if _Hpols[0].stop - _Hpols[0].start == 1: self._Hpol = _Hpols[1] else: self._Hpol = _Hpols[2] return x, y, z
[docs] def update_E(self): """Add the source to the electric field""" q = self.grid.time_steps_passed vect = self.profile * sin(2 * pi * q / self.period + self.phase_shift) self.grid.E[self.x, self.y, self.z, self._Epol] = vect
[docs] def update_H(self): """Add the source to the magnetic field""" q = self.grid.time_steps_passed vect = self.profile * sin(2 * pi * q / self.period + self.phase_shift) self.grid.H[self.x, self.y, self.z, self._Hpol] = vect
def __repr__(self): return ( f"{self.__class__.__name__}(period={self.period}, " f"amplitude={self.amplitude}, phase_shift={self.phase_shift}, " f"name={repr(self.name)}, polarization={repr(self.polarization)})" ) def __str__(self): s = " " + repr(self) + "\n" x = f"[{self.x.start}, ... , {self.x.stop}]" y = f"[{self.y.start}, ... , {self.y.stop}]" z = f"[{self.z.start}, ... , {self.z.stop}]" s += f" @ x={x}, y={y}, z={z}\n" return s
[docs]class SoftArbitraryPointSource: r""" A source placed at a single point (grid cell) in the grid. This source is special: it's both a source and a detector. Unlike the other sources, the input is a voltage, not an electric field. (really? why? should we convert back and forth?) For electrical measurements I've only needed a single-index source, so I don't know how the volume/line sources above work. We want the FFT function to operate over any detector. Maybe all sources should take an arbitary waveform argument? Each index in the *waveform* array represents 1 value at a timestep. There are many different *geometries* of "equivalent sources". The detector/source paradigm used in /fdtd might perhaps not correspond to this in an ideal fashion. It's not intuitively clear to me what a "soft" source would imply in the optical case, or what impedance even means for a laser. /fdtd/ seems to have found primary use in optical circles, so the default Z should probably be 0. "Whilst established for microwaves and electrical circuits, this concept has only very recently been observed in the optical domain, yet is not well defined or understood."[1] [1]: Optical impedance of metallic nano-structures, M. Mazilu and K. Dholakia https://doi.org/10.1364/OE.14.007709 [2]: http://www.gwoptics.org/learn/02_Plane_waves/01_Fabry_Perot_cavity/02_Impedance_matched.php /\/\- """
[docs] def __init__( self, waveform_array: ndarray, name: str = None, impedance: float = 0.0 ): """Create Args: waveform_array """ self.grid = None self.name = name self.current_detector = None self.waveform_array = waveform_array self.impedance = impedance self.input_voltage = [] # voltage hard-imposed by the source self.source_voltage = [] #
# "field" rather than "voltage" might be more meaningful # FIXME: these voltage time histories have a different dimensionality def _register_grid(self, grid: Grid, x: Number, y: Number, z: Number): """Register a grid for the source. Args: grid: the grid to place the source into. x: The x-location of the source in the grid y: The y-location of the source in the grid z: The z-location of the source in the grid Note: As its name suggests, this source is a POINT source. Hence it should be placed at a single coordinate tuple int the grid. """ self.grid = grid self.grid.sources.append(self) if self.name is not None: if not hasattr(grid, self.name): setattr(grid, self.name, self) else: raise ValueError( f"The grid already has an attribute with name {self.name}" ) try: (x,), (y,), (z,) = x, y, z except (TypeError, ValueError): raise ValueError("a point source should be placed on a single grid cell.") self.x, self.y, self.z = grid._handle_tuple((x, y, z)) if self.name is not None: detector_name += "_I" else: detector_name = None self.current_detector = CurrentDetector(name=detector_name) grid[x, y, z] = self.current_detector
[docs] def update_E(self): # It is important that this step happen between the E-field update and the # H-field update. if self.grid.time_steps_passed < self.waveform_array.shape[0]: # check for off-by-one error here input_voltage = self.waveform_array[self.grid.time_steps_passed] else: input_voltage = 0.0 # one could taper the last value off smoothly instead if self.grid.time_steps_passed > 0: current = self.current_detector.I[-1][0][0][0] else: current = 0.0 if self.impedance > 0: source_resistive_voltage = self.impedance * current output_voltage = input_voltage + source_resistive_voltage else: output_voltage = input_voltage # right now, this does not compensate for the cell's permittivity! self.grid.E[self.x, self.y, self.z, 2] += ( output_voltage / self.grid.grid_spacing ) self.input_voltage.append([[[input_voltage]]]) self.source_voltage.append([[[output_voltage]]])
[docs] def update_H(self): pass
def __repr__(self): return f"{self.__class__.__name__}()" def __str__(self): s = " " + repr(self) + "\n" x = f"{self.x}" y = f"{self.y}" z = f"{self.z}" s += f" @ x={x}, y={y}, z={z}\n" return s