flaport / fdtd

A 3D electromagnetic FDTD simulator written in Python with optional GPU support
https://fdtd.readthedocs.io
MIT License
454 stars 116 forks source link

Pulsing field of simple point sources? #73

Open ss32 opened 1 month ago

ss32 commented 1 month ago

I'm running a modified version of the first example and am wondering why the whole field pulses at what I assume is the sim unit period of the sources. Shouldn't the sources pulse at 1 "Hz" and the field's amplitude remain more or less constant (standing wave) once the waves have propagated throughout the simulation space?

test

import matplotlib.pyplot as plt

import fdtd
import fdtd.backend as bd

fdtd.set_backend("torch.cuda")

WAVELENGTH = 1550e-9
C: float = 299_792_458.0
FREQUENCY = C / WAVELENGTH

grid_x_meters = 2.5e-5
grid_y_meters = 1.5e-5
spacing = 0.1 * WAVELENGTH
grid_x_pixels = grid_x_meters / spacing
grid_y_pixels = grid_y_meters / spacing
# Create ftdt grid
grid = fdtd.Grid(
    (grid_x_meters, grid_y_meters, 1),
    grid_spacing=spacing,
    permittivity=1.0,
    permeability=1.0,
)

source_x = 8.525e-6
source_y = 1.1625e-05
source_z = 0

print(f"Frequency: {FREQUENCY:.2e} Hz")
# Boundary conditions
# x bounds
grid[0:10, :, :] = fdtd.PML(name="pml_xlow")
grid[-10:, :, :] = fdtd.PML(name="pml_xhigh")
# y bounds
grid[:, 0:10, :] = fdtd.PML(name="pml_ylow")
grid[:, -10:, :] = fdtd.PML(name="pml_yhigh")

grid[:, :, -1] = fdtd.PeriodicBoundary(name="zbounds")

# Add sources
grid[source_x, source_y, source_z] = fdtd.PointSource(period=WAVELENGTH / C, name="pointsource1")
grid[source_x, source_y - 8e-6, source_z] = fdtd.PointSource(
    period=WAVELENGTH / C ,
    name="pointsource2",
)

# Add detector
grid[source_x-1e-6, :, 0] = fdtd.LineDetector(name="detector")

# Add an object
grid[11:32, 30:84, 0:1] = fdtd.AnisotropicObject(permittivity=2.5, name="object")

# Run the sim
# grid.run(total_time=500, progress_bar=True)
for i in range(500):
    grid.step()
    grid.visualize(z=0, animate=True, index=i, save=True, folder="test")
    plt.title(f"{i:3.0f}")
    plt.show()