Poisson FDM Solver

A simple framework for setting up electrostatic problems and solving them using the finite difference method.

Basics

In order to set up an electrostatic problem, the grid values of ρ, ε and the boundary conditions can be set.

drawing_ops contains functions for drawing basic geometric objects onto a grid. When an input argument is passed to the function, the specified object is drawn onto it.

import drawing_ops as ops
N = 256
U = np.zeros((N, N))
U = ops.plate_capacitor(U, center=(N/2, N/2),
                           length=N/2,
                           distance=N/8,
                           rotation=90)

Alternatively, if only geometric properties are specified, the function returns another function that can draw the object onto an input.

N = 256
capacitor = ops.plate_capacitor(center=(N/2, N/2), 
                                length=N/2, 
                                distance=N/8, 
                                rotation=90)
U = capacitor(np.zeros((N, N)))

The function solve_poisson receives the grid values as input and returns the finite-difference solution of φ.

phi = solve_poisson(U)

(1) Plate capacitor

capacitor = ops.plate_capacitor(center=(N/2, N/2), 
                                length=N/2, 
                                distance=N/8, 
                                rotation=90)
U = capacitor(np.zeros((N, N)))
phi = solve_poisson(U)
plot_images([U, phi], ['U', 'φ'], cmap='afmhot')

Ey, Ex = np.gradient(phi)
Ex, Ey = -Ex, -Ey
E = np.sqrt(Ex**2 + Ey**2)

(2) Dielectric block between two plates

capacitor = ops.plate_capacitor(center=(N/2-1, N/2-1), 
                                length=N/1.5, 
                                distance=N/8, 
                                rotation=90, 
                                values=[-10, 10])
U = capacitor(np.zeros((N, N)))

Eps = np.ones((N, N))
Eps = ops.rectangle(Eps, center=(N/2-1, N/2-1), 
                         wh=(N//16, N//16), 
                         filled=True, 
                         value=5)

phi = solve_poisson(U=U, Eps=Eps)

(3) Two oppositely charged circles

charges = ops.composition(
    ops.circle(center=(N/2-N/5, N/2), radius=N/16, filled=True, value=1),
    ops.circle(center=(N/2+N/5, N/2), radius=N/16, filled=True, value=-1)
)
Rho = charges(np.zeros((N, N)))

phi = solve_poisson(U=np.zeros((N, N)), Rho=Rho)

(4) Charged circle in front of a plate

Rho = ops.circle(np.zeros((N, N)), center=(N/2+N/8, N/2), 
                                   radius=N/32, 
                                   filled=True)
U = ops.line(np.zeros((N, N)), center=(N/2-N/8, N/2), 
                               length=N/2, 
                               rotation=90, 
                               value=1)

phi = solve_poisson(U=U, Rho=Rho)

(5) Charged circle in front of a dielectric half-space

Eps = np.ones((N, N))
Eps[:,:N//2-N//8] = 3
Rho = ops.circle(np.zeros((N, N)), center=(N/2+N/16, N/2), 
                                   radius=N/32, 
                                   filled=True)

phi = solve_poisson(Eps=Eps, Rho=Rho)

Make videos

frames = make_frames(ts, U, stars)
frames_c = colorize_frames(frames['E'], num_levels=30)

plt.figure(figsize=(10, 10))
plot_frames(frames_c, num=10, nrow=5)

make_video('data/E_stars_512_fps20_l30_m15.mp4', frames_c, fps=20, fourcc='h264')

Examples

distance = lambda t: N/5 + N/10*np.sin(2*np.pi*t)
cap = lambda t: ops.plate_capacitor(center=(N/2, N/2), 
                                    length=N/2, 
                                    distance=distance(t), 
                                    values=[-10, 10], 
                                    plate_width=N//50,
                                    rotation=90)

E_cap_512_fps20_lN_m15.mp4

theta = lambda t: 360*t
arcs2 = lambda t: ops.composition(
  ops.arc(center=(N/2, N/2), radius=N/3, start=0, end=250, rotation=theta(t), value=10),
  ops.arc(center=(N/2, N/2), radius=N/3-N/10, start=0, end=250, rotation=45-theta(t), value=-10),
  ops.arc(center=(N/2, N/2), radius=N/3-2*N/10, start=0, end=250, rotation=90+theta(t), value=10),
  ops.arc(center=(N/2, N/2), radius=N/3-3*N/10, start=0, end=250, rotation=135-theta(t), value=-10)
)

E_arcs2_512_20_N.mp4

voltage = lambda t, phase: 10*np.sin(2*np.pi*t + phase/180*np.pi)                                        
stars = lambda t: ops.composition(
  ops.star(center=(N/2-N/5, N/2+N/5), N_angles=6, ro=N/10, value=voltage(t, phase=0)),                   
  ops.star(center=(N/2+N/5, N/2+N/5), N_angles=6, ro=N/10, value=voltage(t, phase=120)),                        
  ops.star(center=(N/2, N/2-N/5), N_angles=6, ro=N/10, value=voltage(t, phase=-120)))

E_stars_512_fps20_l30_m15.mp4

References

Nagel 2012 - Solving the Generalized Poisson Equation Using the Finite-Difference Method (FDM)