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

moving source #30

Closed GenB31415 closed 2 years ago

GenB31415 commented 2 years ago

Dear author, many thanks for your interesting package! In my problem I have to study the field of a source moving with the uniform velocity along the grid. Is there any way to indicate the current position of the point source (depending on current time) while the field update? Thank you in advance.

flaport commented 2 years ago

Hey @GenB31415 ,

Moving sources is not really supported by this package. However you can probably hack it as follows:


# define the position x, y and z of the source at each time t here:
x = ...
y = ...
z = ...
t = ...

for i, t in enumerate(t):
    grid.step()
    source = grid.sources.pop(-1) # make sure the source is no longer part of the grid
    grid[x[i], y[i], z[i]] = source
GenB31415 commented 2 years ago

Dear author, many thanks for the idea.I did that but I the error is generated: The grid already has an attribute with name pointsource My code is bellow (slightly modified code 01-basic-example.py ). Please let me know how to overcome this issue. Thank you. import matplotlib.pyplot as plt import fdtdimport fdtd.backend as bdfdtd.set_backend("numpy")# ConstantsWAVELENGTH = 1550e-9SPEED_LIGHT: float = 299_792_458.0  # [m/s] speed of light# Simulation# create FDTD Gridgrid = fdtd.Grid(    (2.5e-5, 1.5e-5, 1),    grid_spacing=0.1 * WAVELENGTH,    permittivity=1.0,    permeability=1.0,)

boundaries# grid[0, :, :] = fdtd.PeriodicBoundary(name="xbounds")grid[0:10, :, :] = fdtd.PML(name="pml_xlow")grid[-10:, :, :] = fdtd.PML(name="pml_xhigh")

grid[:, 0, :] = fdtd.PeriodicBoundary(name="ybounds")grid[:, 0:10, :] = fdtd.PML(name="pml_ylow")grid[:, -10:, :] = fdtd.PML(name="pml_yhigh")

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

sources# grid[50:55, 70:75, 0] = fdtd.LineSource(#     period=WAVELENGTH / SPEED_LIGHT, name="linesource"# )#grid[100, 60, 0] = fdtd.PointSource(grid[120, 50, 0] = fdtd.PointSource(    period=WAVELENGTH / SPEED_LIGHT, name="pointsource",)# detectors

grid[12e-6, :, 0] = fdtd.LineDetector(name="detector")# objects

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

*- - - - - - start - - - - -** - this part is added !- -# define the position x, y and z of the source at each time t here:x0 = 120y0 = 50z0 = 0tFin=50t = range(0,tFin)velocity = 0.1 # <---- velocity of a point source only in x direction

for i, t in enumerate(t):    grid.step()        posX = x0-velocity*t # <---- move source in next position x at time increase    posY = y0 # <---- fixed position y    posZ = z0 # <---- fixed position t    source = grid.sources.pop(-1) # make sure the source is no longer part of the grid    grid[posX, posY, posZ] = source # <--- Error: The grid already has an attribute with name pointsource    # *** - - - - - end- - - - - - - - -

Run simulationgrid.run(tFin, progress_bar=False)# Visualizationfig, axes = plt.subplots(2, 3, squeeze=False)titles = ["Ex: xy", "Ey: xy", "Ez: xy", "Hx: xy", "Hy: xy", "Hz: xy"]

fields = bd.stack(    [        grid.E[:, :, 0, 0],        grid.E[:, :, 0, 1],        grid.E[:, :, 0, 2],        grid.H[:, :, 0, 0],        grid.H[:, :, 0, 1],        grid.H[:, :, 0, 2],    ]) m = max(abs(fields.min().item()), abs(fields.max().item())) for ax, field, title in zip(axes.ravel(), fields, titles):    ax.set_axis_off()    ax.set_title(title)    ax.imshow(bd.numpy(field), vmin=-m, vmax=m, cmap="RdBu") plt.show() plt.figure()grid.visualize(z=0)

On Monday, August 23, 2021, 09:21:48 AM CDT, Floris Laporte ***@***.***> wrote:  

Hey @GenB31415 ,

Moving sources is not really supported by this package. However you can probably hack it as follows:

define the position x, y and z of the source at each time t here:

x = ... y = ... z = ... t = ...

for i, t in enumerate(t): grid.step() source = grid.sources.pop(-1) # make sure the source is no longer part of the grid grid[x[i], y[i], z[i]] = source — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android.

GenB31415 commented 2 years ago

(repetition) Dear @flaport , many thanks for the idea. I did that but I the error is generated: The grid already has an attribute with name pointsource

My code is bellow (slightly modified code 01-basic-example.py ). Please let me know how to overcome this issue. Thank you.

#fdtd.readthedocs.io/en/latest/examples/01-basic-example.html
import matplotlib.pyplot as plt
import fdtd
import fdtd.backend as bd
fdtd.set_backend("numpy")
# Constants
WAVELENGTH = 1550e-9
SPEED_LIGHT: float = 299_792_458.0  # [m/s] speed of light
# Simulation
# create FDTD Grid
grid = fdtd.Grid(
    (2.5e-5, 1.5e-5, 1),
    grid_spacing=0.1 * WAVELENGTH,
    permittivity=1.0,
    permeability=1.0,
)
# boundaries
# grid[0, :, :] = fdtd.PeriodicBoundary(name="xbounds")
grid[0:10, :, :] = fdtd.PML(name="pml_xlow")
grid[-10:, :, :] = fdtd.PML(name="pml_xhigh")

# grid[:, 0, :] = fdtd.PeriodicBoundary(name="ybounds")
grid[:, 0:10, :] = fdtd.PML(name="pml_ylow")
grid[:, -10:, :] = fdtd.PML(name="pml_yhigh")

grid[:, :, 0] = fdtd.PeriodicBoundary(name="zbounds")
# sources
# grid[50:55, 70:75, 0] = fdtd.LineSource(
#     period=WAVELENGTH / SPEED_LIGHT, name="linesource"
# )
#grid[100, 60, 0] = fdtd.PointSource(
grid[120, 50, 0] = fdtd.PointSource(
    period=WAVELENGTH / SPEED_LIGHT, name="pointsource",
)
# detectors
#grid[12e-6, :, 0] = fdtd.LineDetector(name="detector")
# objects
grid[11:32, 30:84, 0:1] = fdtd.AnisotropicObject(permittivity=2.5, name="object")
# - - - - - - start - - - - - - - -
# define the position x, y and z of the source at each time t here:
x0 = 120
y0 = 50
z0 = 0
tFin=50
t = range(0,tFin)
velocity = 0.1 # <---- velocity of a point source only in x direction

for i, t in enumerate(t):
    grid.step()    
    posX = x0-velocity*t # <---- move source in next position x at time increase
    posY = y0 # <---- fixed position y
    posZ = z0 # <---- fixed position t
    source = grid.sources.pop(-1) # make sure the source is no longer part of the grid
    grid[posX, posY, posZ] = source    
# - - - - - end- - - - - - - - -

# Run simulation
grid.run(tFin, progress_bar=False)
# Visualization
fig, axes = plt.subplots(2, 3, squeeze=False)
titles = ["Ex: xy", "Ey: xy", "Ez: xy", "Hx: xy", "Hy: xy", "Hz: xy"]

fields = bd.stack(
    [
        grid.E[:, :, 0, 0],
        grid.E[:, :, 0, 1],
        grid.E[:, :, 0, 2],
        grid.H[:, :, 0, 0],
        grid.H[:, :, 0, 1],
        grid.H[:, :, 0, 2],
    ]
)
m = max(abs(fields.min().item()), abs(fields.max().item()))

for ax, field, title in zip(axes.ravel(), fields, titles):
    ax.set_axis_off()
    ax.set_title(title)
    ax.imshow(bd.numpy(field), vmin=-m, vmax=m, cmap="RdBu")

plt.show()
plt.figure()
grid.visualize(z=0)
flaport commented 2 years ago

You gave your source the name 'pointsource', so it's available as an attribute to the grid, i.e. grid.pointsource. In that case you also have to delete that attribute in the loop:

source = grid.sources.pop(-1)
del grid.pointsource
GenB31415 commented 2 years ago

Hey @flaport, thanks, now it works. The code (run or uncomment corresponding line)

import matplotlib.pyplot as plt
import fdtd
import fdtd.backend as bd
fdtd.set_backend("numpy")

OneMicron=1e-6
WAVELENGTH = OneMicron 
LengthScale = OneMicron
SPEED_LIGHT: float = 299_792_458.0  # [m/s] speed of light
indRefr=1.15
grid = fdtd.Grid(
    (25*LengthScale, 25*LengthScale, 1),  
    grid_spacing=0.1 * LengthScale, permittivity=indRefr**2, permeability=1.0,
)
grid[0:10, :, :] = fdtd.PML(name="pml_xlow")
grid[-10:, :, :] = fdtd.PML(name="pml_xhigh")    
grid[:, 0:10, :] = fdtd.PML(name="pml_ylow")
grid[:, -10:, :] = fdtd.PML(name="pml_yhigh")
x0 = 17.0*LengthScale; y0 = 12.5*LengthScale; z0 = 0; 
periodTime= LengthScale / SPEED_LIGHT
grid[x0, y0, z0] = fdtd.PointSource( period = periodTime )

# velocity of a point source only in x direction
#velocity = 0.69; tFin=20*1  
#velocity = 0.6; tFin=20*1  
#velocity = 0.5; tFin=20*1  
#velocity = 0.3; tFin=20*2  
velocity = 0.2; tFin=20*3  
#velocity = 0.1; tFin=20*4
#velocity = 0.05; tFin=20*4
#velocity = 0.005; tFin=20*5

tAll = range(0, tFin)  
for i, t in enumerate(tAll):
    grid.step()    
    tt = t*LengthScale      
    posX = (x0-velocity*tt) # <---- move source in next position x when time increases
    posY = y0               # <---- fixed position y
    posZ = z0               # <---- fixed position z
    source = grid.sources.pop(-1) # make sure the source is no longer part of the grid    
    source = fdtd.PointSource(period=periodTime, name=(f"Cher: i={i}"))
    grid[posX, posY, posZ] = source    

grid.run(tFin, progress_bar=False)
fig, axes = plt.subplots(2, 3, squeeze=False)
titles = ["Ex: xy", "Ey: xy", "Ez: xy", "Hx: xy", "Hy: xy", "Hz: xy"]
fields = bd.stack(
    [
        grid.E[:, :, 0, 0],
        grid.E[:, :, 0, 1],
        grid.E[:, :, 0, 2],
        grid.H[:, :, 0, 0],
        grid.H[:, :, 0, 1],
        grid.H[:, :, 0, 2],
    ]
)
m = max(abs(fields.min().item()), abs(fields.max().item()))
import numpy as np
for ax, field, title in zip(axes.ravel(), fields, titles):
    ax.set_axis_off()
    ax.set_title(title)    
    ax.imshow(bd.numpy(field), vmin=-m, vmax=m, cmap="hsv")    

plt.show()
plt.figure()
plt.title(source.name +f", v={velocity}, n={indRefr}")
grid.visualize(z=0,cmap="nipy_spectral_r",)
flaport commented 2 years ago

great to hear it's working now :)