Maxwell solver in two dimensions with FDTD scheme

Contents

Maxwell solver in two dimensions with FDTD scheme#

\[ \frac{\partial H_z}{\partial t} = \frac{\partial E_x}{\partial y} - \frac{\partial E_y}{\partial x} ;\qquad \frac{\partial E_x}{\partial t} = \frac{\partial H_z}{\partial y} ;\qquad \frac{\partial E_y}{\partial t} = - \frac{\partial H_z}{\partial x} \]

fdtd $$ H_z \big|^{n+1/2}{i+1/2,j+1/2} = H_z \big|^{n-1/2}{i+1/2,j+1/2} + \frac{dt}{dy} \big(E_x \big|^{n}{i+1/2,j+1} - E_x \big|^{n}{i+1/2,j} \big)

  • \frac{dt}{dx} \big( E_y \big|^{n}{i+1,j+1/2} - E_y \big|^{n}{i,j+1/2} \big) $$

\[ E_x \big|^{n+1}_{i+1/2,j} = E_x \big|^{n}_{i+1/2,j} + \frac{dt}{dy} \big( H_z \big|^{n+1/2}_{i+1/2,j+1/2} - H_z \big|^{n+1/2}_{i-1/2, j-1/2} \big) \]
\[ E_y \big|^{n+1}_{i,j+1/2} = E_y \big|^{n}_{i,j+1/2} - \frac{dt}{dx} \big( H_z \big|^{n+1/2}_{i+1/2,j+1/2} - H_z \big|^{n+1/2}_{i-1/2, j+1/2} \big) \]

Description of the scheme

%matplotlib inline
%config InlineBackend.figure_format = 'retina'

import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import axes3d
import matplotlib.animation as animation
from IPython.display import HTML

plt.rcParams['figure.figsize'] = (10,6)
# Mesh parameters
nx, ny = 101, 101
vx, dx = np.linspace(0, 1, nx, endpoint=True, retstep=True)
vy, dy = np.linspace(0, 1, ny, endpoint=True, retstep=True)

#Initialize Ex, Ey when time = 0
ex = np.zeros((nx-1, ny), dtype=np.double)  
ey = np.zeros((nx, ny-1), dtype=np.double) 
nbiter = 500   # time loop size
dt = 0.001     # time step
m, n = 2, 2
omega = np.sqrt((m*np.pi)**2+(n*np.pi)**2)
# Create the staggered grid for Bz
x, y = np.meshgrid(0.5*(vx[:-1]+vx[1:]), 0.5*(vy[:-1]+vy[1:]))
fig = plt.figure()
ax = axes3d.Axes3D(fig)

#Initialize Bz when time = - dt / 2
hz = - np.cos(m*np.pi*y) * np.cos(n*np.pi*x) * np.cos(omega*(-0.5*dt))
wframe = ax.plot_wireframe(x, y, hz, rstride=2, cstride=2)
ax.set_zlim(-1,1);
<Figure size 1000x600 with 0 Axes>

numpy#

def faraday( ex, ey, hz ) : 
    "faraday equation Bz(t+dt/2) -> Bz(t-dt/2) + dt f(E(t))"
    return hz + dt * ((ex[:, 1:]-ex[:, :-1]) / dy - (ey[1:, :]-ey[:-1, :]) / dx)

def ampere_maxwell( hz, ex, ey):
    " Ampere-Maxwell equation E(t+dt) -> E(t) + dt g(Bz(t+dt/2)) "
    ex[:, 1:-1] += dt*(hz[:, 1:]-hz[:, :-1]) / dy
    ey[1:-1, :] += - dt*(hz[1:, :]-hz[:-1, :]) / dx

    # periodic boundary conditions
    ex[:, 0] += dt*(hz[:, 0]-hz[:, -1]) / dy
    ex[:, -1] = ex[:, 0]
    ey[0, :] += - dt*(hz[0, :]-hz[-1, :]) / dx
    ey[-1, :] = ey[0, :]
    
    return ex, ey
def update(i, ax, fig):
    ax.cla()

    global ex, ey, hz

    hz = faraday( ex, ey, hz)
    ex, ey = ampere_maxwell( hz, ex, ey)
    
    wframe = ax.plot_wireframe(x, y, hz, rstride=2, cstride=2)
    ax.set_zlim(-1, 1)
    return wframe,
ani = animation.FuncAnimation(fig, update,
                              frames=range(200),
                              fargs=(ax, fig), interval=100)
%%time
HTML(ani.to_html5_video())
CPU times: user 2.39 s, sys: 102 ms, total: 2.5 s
Wall time: 3.27 s
%%time

from tqdm.notebook import tqdm

nx, ny = 512, 512
vx, dx = np.linspace(0, 1, nx, endpoint=True, retstep=True)
vy, dy = np.linspace(0, 1, ny, endpoint=True, retstep=True)

ex = np.zeros((nx-1, ny), dtype=np.double)  
ey = np.zeros((nx, ny-1), dtype=np.double) 
dt = 0.001     # time step
m, n = 2, 2
omega = np.sqrt((m*np.pi)**2+(n*np.pi)**2)
x, y = np.meshgrid(0.5*(vx[:-1]+vx[1:]), 0.5*(vy[:-1]+vy[1:]))

hz = - np.cos(m*np.pi*y) * np.cos(n*np.pi*x) * np.cos(omega*(-0.5*dt))

for t in tqdm(range(1000)):
    
    hz = faraday( ex, ey, hz)
    ex, ey = ampere_maxwell( hz, ex, ey)
    
CPU times: user 2.37 s, sys: 3.83 s, total: 6.2 s
Wall time: 6.18 s
%load_ext fortranmagic

fortran#

%%fortran 

subroutine faraday_fortran( ex, ey, bz, dx, dy, dt, nx, ny)
implicit none

real(8), intent(in) :: ex(nx-1,ny)
real(8), intent(in) :: ey(nx,ny-1)
real(8), intent(inout) :: bz(nx-1,ny-1)
integer, intent(in) :: nx, ny
real(8), intent(in) :: dx, dy, dt

integer :: i, j
real(8) :: dex_dx, dey_dy
real(8) :: dex_dy, dey_dx
        
do j=1,ny-1
do i=1,nx-1
   dex_dy  = (ex(i,j+1)-ex(i,j)) / dy
   dey_dx  = (ey(i+1,j)-ey(i,j)) / dx
   bz(i,j) = bz(i,j) +  dt * (dex_dy - dey_dx)
end do
end do

end subroutine faraday_fortran
Traceback (most recent call last):
  File "<frozen runpy>", line 198, in _run_module_as_main
  File "<frozen runpy>", line 88, in _run_code
  File "/usr/share/miniconda3/envs/runenv/lib/python3.13/site-packages/numpy/f2py/__main__.py", line 5, in <module>
    main()
    ~~~~^^
  File "/usr/share/miniconda3/envs/runenv/lib/python3.13/site-packages/numpy/f2py/f2py2e.py", line 781, in main
    run_compile()
    ~~~~~~~~~~~^^
  File "/usr/share/miniconda3/envs/runenv/lib/python3.13/site-packages/numpy/f2py/f2py2e.py", line 753, in run_compile
    builder.compile()
    ~~~~~~~~~~~~~~~^^
  File "/usr/share/miniconda3/envs/runenv/lib/python3.13/site-packages/numpy/f2py/_backends/_meson.py", line 192, in compile
    self.run_meson(self.build_dir)
    ~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^
  File "/usr/share/miniconda3/envs/runenv/lib/python3.13/site-packages/numpy/f2py/_backends/_meson.py", line 185, in run_meson
    self._run_subprocess_command(setup_command, build_dir)
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/share/miniconda3/envs/runenv/lib/python3.13/site-packages/numpy/f2py/_backends/_meson.py", line 181, in _run_subprocess_command
    subprocess.run(command, cwd=cwd, check=True)
    ~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/share/miniconda3/envs/runenv/lib/python3.13/subprocess.py", line 554, in run
    with Popen(*popenargs, **kwargs) as process:
         ~~~~~^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/share/miniconda3/envs/runenv/lib/python3.13/subprocess.py", line 1039, in __init__
    self._execute_child(args, executable, preexec_fn, close_fds,
    ~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
                        pass_fds, cwd, env,
                        ^^^^^^^^^^^^^^^^^^^
    ...<5 lines>...
                        gid, gids, uid, umask,
                        ^^^^^^^^^^^^^^^^^^^^^^
                        start_new_session, process_group)
                        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/share/miniconda3/envs/runenv/lib/python3.13/subprocess.py", line 1972, in _execute_child
    raise child_exception_type(errno_num, err_msg, err_filename)
FileNotFoundError: [Errno 2] No such file or directory: 'meson'
Cannot use distutils backend with Python>=3.12, using meson backend instead.
Using meson backend
Will pass --lower to f2py
See https://numpy.org/doc/stable/f2py/buildtools/meson.html
Reading fortran codes...
	Reading file '/home/runner/.cache/ipython/fortranmagic/d694754d/_fortran_magic_348c2b611ed5c55864f30bdb3f0d267a.f90' (format:free)
Post-processing...
	Block: _fortran_magic_348c2b611ed5c55864f30bdb3f0d267a
			Block: faraday_fortran
Applying post-processing hooks...
  character_backward_compatibility_hook
Post-processing (stage 2)...
Building modules...
    Building module "_fortran_magic_348c2b611ed5c55864f30bdb3f0d267a"...
    Generating possibly empty wrappers"
    Maybe empty "_fortran_magic_348c2b611ed5c55864f30bdb3f0d267a-f2pywrappers.f"
        Constructing wrapper function "faraday_fortran"...
          faraday_fortran(ex,ey,bz,dx,dy,dt,[nx,ny])
    Wrote C/API module "_fortran_magic_348c2b611ed5c55864f30bdb3f0d267a" to file "./_fortran_magic_348c2b611ed5c55864f30bdb3f0d267amodule.c"
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[10], line 1
----> 1 get_ipython().run_cell_magic('fortran', '', '\nsubroutine faraday_fortran( ex, ey, bz, dx, dy, dt, nx, ny)\nimplicit none\n\nreal(8), intent(in) :: ex(nx-1,ny)\nreal(8), intent(in) :: ey(nx,ny-1)\nreal(8), intent(inout) :: bz(nx-1,ny-1)\ninteger, intent(in) :: nx, ny\nreal(8), intent(in) :: dx, dy, dt\n\ninteger :: i, j\nreal(8) :: dex_dx, dey_dy\nreal(8) :: dex_dy, dey_dx\n\ndo j=1,ny-1\ndo i=1,nx-1\n   dex_dy  = (ex(i,j+1)-ex(i,j)) / dy\n   dey_dx  = (ey(i+1,j)-ey(i,j)) / dx\n   bz(i,j) = bz(i,j) +  dt * (dex_dy - dey_dx)\nend do\nend do\n\nend subroutine faraday_fortran\n')

File /usr/share/miniconda3/envs/runenv/lib/python3.13/site-packages/IPython/core/interactiveshell.py:2549, in InteractiveShell.run_cell_magic(self, magic_name, line, cell)
   2547 with self.builtin_trap:
   2548     args = (magic_arg_s, cell)
-> 2549     result = fn(*args, **kwargs)
   2551 # The code below prevents the output from being displayed
   2552 # when using magics with decorator @output_can_be_silenced
   2553 # when the last Python token in the expression is a ';'.
   2554 if getattr(fn, magic.MAGIC_OUTPUT_CAN_BE_SILENCED, False):

File /usr/share/miniconda3/envs/runenv/lib/python3.13/site-packages/fortranmagic.py:516, in FortranMagics.fortran(self, line, cell)
    511 res = self._run_f2py(f2py_args + ['-m', module_name, '-c',
    512                                   f_f90_file],
    513                      verbosity=args.verbosity,
    514                      fflags=fflags)
    515 if res != 0:
--> 516     raise RuntimeError("f2py failed, see output")
    518 self._code_cache[key] = module_name
    519 module = _imp_load_dynamic(module_name, module_path)

RuntimeError: f2py failed, see output
%%fortran

subroutine amperemaxwell_fortran(ex, ey, bz, dx, dy, dt, nx, ny)

    implicit none
    integer, intent(in):: nx, ny
    real(8), intent(in):: dx, dy, dt
    real(8), dimension(nx-1, ny-1), intent(inout) :: bz
    real(8), dimension(nx-1, ny),   intent(inout) :: ex
    real(8), dimension(nx, ny-1),   intent(inout) :: ey
    integer:: i, j
    real(8):: dbz_dx, dbz_dy
    real(8), parameter:: csq = 1d0
    
    do i = 1, nx-1
       dbz_dy = (bz(i, 1)-bz(i, ny-1)) / dy ! periodic BC
       ex(i, 1) = ex(i, 1) + dt*csq*dbz_dy
       ex(i, ny) = ex(i, 1)
    end do
    
    do j = 1, ny-1
       dbz_dx = (bz(1,j)-bz(nx-1,j)) / dx  ! periodic BC
       ey(1,j) = ey(1,j) - dt*csq*dbz_dx
       ey(nx,j) = ey(1,j)
    end do
    
    do j=2,ny-1
        do i=1,nx-1
           dbz_dy = (bz(i,j)-bz(i,j-1)) / dy
           ex(i,j) = ex(i,j) + dt*csq*dbz_dy 
        end do
    end do
    
    do j=1,ny-1
        do i=2,nx-1
           dbz_dx = (bz(i,j)-bz(i-1,j)) / dx
           ey(i,j) = ey(i,j) - dt*csq*dbz_dx 
        end do
    end do

end subroutine amperemaxwell_fortran
%%time

from tqdm.notebook import tqdm

ex.fill(0.0)
ey.fill(0.0)
hz = - np.cos(m*np.pi*y) * np.cos(n*np.pi*x) * np.cos(omega*(-0.5*dt))
ex = np.asfortranarray(ex)
ey = np.asfortranarray(ey)
hz = np.asfortranarray(hz)

for t in tqdm(range(1000)):
    
    faraday_fortran( ex, ey, hz, dx, dy, dt, nx, ny)
    amperemaxwell_fortran(ex, ey, hz, dx, dy, dt, nx, ny)