WI add tests

nidn/fdtd/detectors.py

@@ -14,7 +14,7 @@
 # relative
 from .grid import Grid
 from .backend import backend as bd
-from .constants import X, Y, Z
+from .constants import FDTD_X, FDTD_Y, FDTD_Z
 
 ## Detector
 class LineDetector:
@@ -441,20 +441,22 @@ def single_point_current(self, px, py, pz):
  # option for unit factor here? Units will get very complicated otherwise
 
  current_vector_1 = (
- (self.grid.H[px, py - 1, pz, X]) - (self.grid.H[px, py, pz, X])
+ (self.grid.H[px, py - 1, pz, FDTD_X]) - (self.grid.H[px, py, pz, FDTD_X])
  ) * self.grid.grid_spacing
  current_vector_2 = (
- (self.grid.H[px, py, pz, Y]) - (self.grid.H[px - 1, py, pz, Y])
+ (self.grid.H[px, py, pz, FDTD_Y]) - (self.grid.H[px - 1, py, pz, FDTD_Y])
  ) * self.grid.grid_spacing
 
  current_1 = current_vector_1 + current_vector_2
  # current_1 = float(current.cpu())
 
  current_vector_1 = (
- (self.grid.H[px, py - 1, pz - 1, X]) - (self.grid.H[px, py, pz - 1, X])
+ (self.grid.H[px, py - 1, pz - 1, FDTD_X])
+ - (self.grid.H[px, py, pz - 1, FDTD_X])
  ) * self.grid.grid_spacing
  current_vector_2 += (
- (self.grid.H[px, py, pz - 1, Y]) - (self.grid.H[px - 1, py, pz - 1, Y])
+ (self.grid.H[px, py, pz - 1, FDTD_Y])
+ - (self.grid.H[px - 1, py, pz - 1, FDTD_Y])
  ) * self.grid.grid_spacing
  # current_2 = float(current_2.cpu())
  current_2 = current_vector_1 + current_vector_2

nidn/fdtd_integration/calculate_transmission_reflection_coefficients.py

@@ -1,6 +1,8 @@
 from dotmap import DotMap
 from torch.fft import rfft, rfftfreq
-from torch import sqrt, tensor
+import torch
+from loguru import logger
+
 
 from nidn.fdtd_integration.constants import FDTD_GRID_SCALE
 
@@ -49,6 +51,20 @@ def calculate_transmission_reflection_coefficients(
  reflection_coefficient = _fft(true_reflection, cfg) / _fft(
  reflection_signals[0], cfg
  )
+
+ if transmission_coefficient < 0 or transmission_coefficient > 1:
+ raise ValueError(
+ f"The transmission coefficient is outside of the physical range between 0 and 1"
+ )
+
+ if reflection_coefficient < 0 or reflection_coefficient > 1:
+ raise ValueError(
+ f"The reflection coefficient is outside of the physical range between 0 and 1"
+ )
+ if transmission_coefficient + reflection_coefficient > 1:
+ logger.warning(
+ f"The sum of the transmission and reflection coefficient is greater than 1, which is physically impossible"
+ )
  return transmission_coefficient, reflection_coefficient
 
 
@@ -74,9 +90,11 @@ def _fft(signal, cfg: DotMap):
  tuple[array,array]: fourier frequenices and their corresponding values
  """
  sampling_frequencies = (
- cfg.physical_wavelength_range[0] * FDTD_GRID_SCALE / (sqrt(2) * SPEED_OF_LIGHT)
+ cfg.physical_wavelength_range[0]
+ * FDTD_GRID_SCALE
+ / (torch.sqrt(2) * SPEED_OF_LIGHT)
  )
- tensor_signal = tensor(signal)
+ tensor_signal = torch.tensor(signal)
 
  yf = rfft(tensor_signal)
  xf = rfftfreq(cfg.FDTD_niter, sampling_frequencies)

nidn/fdtd_integration/compute_spectrum_fdtd.py

@@ -1,6 +1,6 @@
 from dotmap import DotMap
 from tqdm import tqdm
-from torch import sqrt, tensor
+import torch
 from loguru import logger
 
 from ..trcwa.get_frequency_points import get_frequency_points
@@ -23,9 +23,10 @@ def compute_spectrum_fdtd(permittivity, cfg: DotMap):
  transmission_spectrum = []
  reflection_spectrum = []
  physical_wavelengths, norm_freq = get_frequency_points(cfg)
- logger.debug("Wavelenghts in spectrum")
+ logger.debug("Wavelenghts in spectrum : ")
  logger.debug(physical_wavelengths)
-
+ logger.debug("Number of layers: ")
+ logger.debug(len(permittivity[0, 0, :, 0]))
  # For each wavelength, calculate transmission and reflection coefficents
  disable_progress_bar = logger._core.min_level >= 20
  for i in tqdm(range(len(physical_wavelengths)), disable=disable_progress_bar):
@@ -59,8 +60,6 @@ def compute_spectrum_fdtd(permittivity, cfg: DotMap):
  )
  transmission_signal.append(transmission_material)
  reflection_signal.append(reflection_material)
- time = [i for i in range(len(transmission_signal[0]))]
-
  # Calculate transmission and reflection coefficients,
  # by using the signals from the free space simulation and the material simulation
  (
@@ -77,7 +76,7 @@ def compute_spectrum_fdtd(permittivity, cfg: DotMap):
  logger.debug("Reflection spectrum")
  logger.debug(reflection_spectrum)
 
- return transmission_spectrum, reflection_spectrum
+ return torch.tensor(transmission_spectrum), torch.tensor(reflection_spectrum)
 
 
 def _get_detector_values(transmission_detector, reflection_detector):
@@ -120,7 +119,7 @@ def _get_abs_value_from_3D_signal(signal):
  abs_value = []
  for i in range(len(signal)):
  abs_value.append(
- sqrt(tensor(signal[i][0] ** 2 + signal[i][1] ** 2 + signal[i][2] ** 2))
+ torch.sqrt(signal[i][0] ** 2 + signal[i][1] ** 2 + signal[i][2] ** 2)
  )
  return abs_value
 
@@ -136,7 +135,7 @@ def _average_along_detector(signal):
  """
  avg = []
  for e in signal:
- s = [tensor(0.0), tensor(0.0), tensor(0.0)]
+ s = torch.zeros(3)
  for p in e:
  s[0] += p[0] / len(e)
  s[1] += p[1] / len(e)

nidn/fdtd_integration/init_fdtd.py

@@ -1,6 +1,15 @@
 from dotmap import DotMap
-import fdtd
 
+from ..fdtd import (
+ AbsorbingObject,
+ set_backend,
+ Grid,
+ PML,
+ PeriodicBoundary,
+ LineSource,
+ LineDetector,
+ PointSource,
+)
 from ..utils.global_constants import EPS_0, SPEED_OF_LIGHT, UNIT_MAGNITUDE
 from .constants import FDTD_GRID_SCALE
 
@@ -17,13 +26,13 @@ def init_fdtd(cfg: DotMap, include_object, wavelength, permittivity):
  Returns:
  fdtd:Grid: Grid with all the added object, ready to be run
  """
- fdtd.set_backend("torch")
+ set_backend("torch")
  scaling = UNIT_MAGNITUDE / (cfg.physical_wavelength_range[0] * FDTD_GRID_SCALE)
  x_grid_size = int(
  scaling * (cfg.FDTD_grid[0] + cfg.N_layers * cfg.PER_LAYER_THICKNESS[0])
  )
  y_grid_size = int(cfg.FDTD_grid[1] * scaling)
- grid = fdtd.Grid(
+ grid = Grid(
  (x_grid_size, y_grid_size, 1),
  grid_spacing=cfg.physical_wavelength_range[0] * FDTD_GRID_SCALE,
  permittivity=1.0,
@@ -38,16 +47,20 @@ def init_fdtd(cfg: DotMap, include_object, wavelength, permittivity):
  cfg.FDTD_reflection_detector_x
  + cfg.N_layers * cfg.PER_LAYER_THICKNESS[0]
  )
+ + cfg.FDTD_min_gridpoints_between_detectors
  ),
  int(cfg.FDTD_reflection_detector_x * scaling),
  )
+ assert cfg.FDTD_pulse_type in ["pulse", "continuous"]
+ use_pulse = cfg.FDTD_pulse_type == "pulse"
+
  grid = _add_source(
  grid,
- int(cfg.FDTD_source[0] * scaling),
- int(cfg.FDTD_source[1] * scaling),
+ int(cfg.FDTD_source_position[0] * scaling),
+ int(cfg.FDTD_source_position[1] * scaling),
  wavelength / SPEED_OF_LIGHT,
- cfg.FDTD_use_pulsesource,
- cfg.FDTD_use_pointsource,
+ use_pulse,
+ cfg.FDTD_source_type,
  )
  if include_object:
  for i in range(cfg.N_layers):
@@ -88,17 +101,17 @@ def _add_boundaries(grid, pml_thickness):
  fdtd.Grid: The grid with the added boundaries
  """
  # Add PML boundary to the left side of the grid
- grid[0:pml_thickness, :, :] = fdtd.PML(name="pml_xlow")
+ grid[0:pml_thickness, :, :] = PML(name="pml_xlow")
  # Add PML boundary to the right side of the grid
- grid[-pml_thickness:, :, :] = fdtd.PML(name="pml_xhigh")
+ grid[-pml_thickness:, :, :] = PML(name="pml_xhigh")
  # Add periodic boundaries at both sides in y-direction
- grid[:, 0, :] = fdtd.PeriodicBoundary(name="ybounds")
+ grid[:, 0, :] = PeriodicBoundary(name="ybounds")
  # Add periodic boundaries on both sides in z-direction. Only applicable for 3D grids
  # grid[:, :, 0] = fdtd.PeriodicBoundary(name="zbounds")
  return grid
 
 
-def _add_source(grid, source_x, source_y, period, use_pulse_source, use_point_source):
+def _add_source(grid, source_x, source_y, period, use_pulse_source, source_type):
  """Add a specified source to the fdtd grid
 
  Args:
@@ -113,17 +126,18 @@ def _add_source(grid, source_x, source_y, period, use_pulse_source, use_point_so
  fdtd.Grid: The grid with the added source
  """
 
- if use_point_source:
- grid[source_x, source_y, 0] = fdtd.PointSource(
+ if source_type == "point":
+ grid[source_x, source_y, 0] = PointSource(
  period=period,
- name="linesource",
+ name="pointsource",
  pulse=use_pulse_source,
  cycle=1,
  hanning_dt=1e-15,
  )
+ elif source_type == "line":
+ grid[source_x, :, 0] = LineSource(period=period, name="linesource")
  else:
- grid[source_x, :, 0] = fdtd.LineSource(period=period, name="linesource")
-
+ raise ValueError(f'FDTD_source_type must be either "line" or "point"')
  return grid
 
 
@@ -139,8 +153,8 @@ def _add_detectors(grid, transmission_detector_x, reflection_detector_x):
  fdtd.Grid: The grid with the added detectors
  """
 
- transmission_detector = fdtd.LineDetector(name="t_detector")
- reflection_detector = fdtd.LineDetector(name="r_detector")
+ transmission_detector = LineDetector(name="t_detector")
+ reflection_detector = LineDetector(name="r_detector")
  grid[transmission_detector_x, :, 0] = transmission_detector
  grid[reflection_detector_x, :, 0] = reflection_detector
  return grid, transmission_detector, reflection_detector
@@ -158,9 +172,9 @@ def _add_object(grid, object_start_x, object_end_x, permittivity, frequency):
  Returns:
  fdtd.Grid: The grid with the added object
  """
- # Not sure whether the conductivity should be relative or absolute, i.e. if it should be multiplied with EPS_0.
+ # Not sure whether the conductivity should be relative or absolute, i.e. if it should be multiplied with EPS_0. Multiplied with 2pi to get w(angular frequency)?
  # Since the permittivity is set to 1 for the free space grid, I'll leave it at an relative value for now. Also, permittivity for object is relative.
- grid[object_start_x:object_end_x, :, :] = fdtd.AbsorbingObject(
+ grid[object_start_x:object_end_x, :, :] = AbsorbingObject(
  permittivity=permittivity.real,
  conductivity=permittivity.imag * frequency,
  )

nidn/plots/plot_spectra.py

@@ -3,7 +3,7 @@
 from matplotlib import pyplot as plt
 from matplotlib.ticker import FormatStrFormatter
 from ..utils.convert_units import freq_to_wl, wl_to_phys_wl
-from ..trcwa.compute_spectrum import compute_spectrum
+from ..utils.compute_spectrum import compute_spectrum
 from ..training.model.model_to_eps_grid import model_to_eps_grid
 
 

nidn/tests/calculate_coefficient_test.py

@@ -0,0 +1,26 @@
+import numpy as np
+from ..fdtd_integration.calculate_transmission_reflection_coefficients import (
+ calculate_transmission_reflection_coefficients,
+)
+from ..utils.load_default_cfg import load_default_cfg
+
+
+def test_calculate_coefficient():
+ cfg = load_default_cfg()
+ time_points = [0.002 * np.pi * i for i in range(1000)]
+ signal_a = 2 * np.sin(time_points)
+ signal_b = np.sin(time_points)
+ signal_array = [signal_a, signal_b]
+ (
+ transmission_coefficient_ms,
+ reflection_coefficient_ms,
+ ) = calculate_transmission_reflection_coefficients(
+ signal_array, signal_array, "mean square", cfg
+ )
+ # TODO: Add test for fft, when the method is complete
+ assert transmission_coefficient_ms - 0.25 == 0
+ assert reflection_coefficient_ms - 0.25 == 0
+
+
+if __name__ == "__main__":
+ test_calculate_coefficient()