%==========================================================================
%
% Solve the 2D Poisson equation using the finite-difference method
% (FDM).
%
% In electromagnetics, the generalized Poisson equation takes on the
% form:
% div ( EPS(r) * grad (V(r)) ) = -RHO(r)
%
% where r is a position vector, RHO is the electrical charge distribution
% and EPS (short for "epsilon") is the electric permittivity. The "div"
% operator is the divergence, while "grad" indicates the gradient.
%
% V = Poisson_FDM_Solver_2D(V,BC,EPS,RHO) returns the 2-D potential
% function V to the initial-valued potential function V_0. Both V and
% V_0 are assumed to have units of VOLTS. The matrix EPS contains the
% dielectric constants (unitlesS) for all points in the domain, while the
% matrix RHO contains the charge density distribution (units of C/m^2).
%
% Note that both EPS and RHO may be either purely REAL (indicating a
% purely static simulation) or COMPLEX numbers (indicating a QUASI-static
% simulation with conductivity). If complex values are used, then the
% output voltage will be complex as well, and must be interpreted as a
% phasor potential.
%
% Also note that in order for the quasi-static assumption to be
% valid, the electric field must be reasonably approximated by
%
% E ~ -grad (V) .
%
% In reality, a time-varying system includes a magnetic vector potential
% that contributes to the E-field. However, as long as this contribution
% is very small, then the simulation results are approximately valid.
%
% BOUNDARY CONDITIONS:
%
% The matrix BC represents the boundary conditions, and specifies
% which points are Dirichlet or Neumann boundaries. All points where BC
% is specified with a 1 are treated as Dirichlet boundaries with fixed
% potentials. This means any inital-valued element of V_0 where the
% corresponding element in BC is a 1 will be treated as a fixed,
% unchanging potential.
% If the point is interior to the simulation domain, then any point
% where BC is not equal to 1 will be solved for by the FD algorithm.
% Conversely, if the point is at the edge of the simulation domain, then
% any value where BC is NOT equal to 1 will be treated as a NEUMANN
% condition. That is, dV/dn = 0, where n is the unit normal to the edge
% of the simulation. This is effectively a repeating boundary, and will
% mimic an exact copy of the simulation domain on the opposite side of
% the boundary.
%
% V = Poisson_FDM_Solver_2D(V,BC,eps,RHO,h) performs the same
% operation, but with the grid spacing h specified in units of meters.
% This only has an effect on the simulation when RHO has nonzero entries.
% Note that the default value is h = 1.0 (meters) if not specified.
%
% GRID SETUP:
%
% It is important to remember that the dielectric samples in EPS are
% staggered from the potentials in V (and RHO) according to the following
% convention:
%
%
% O O O O = Voltage Sample
%
% X X X = Dielectric Sample
%
% O O O
%
% X X
%
% O O O
%
%
% Note how this effectively places the voltage samples at the boundaries
% between dielectric regions. It also means that the size of the EPS
% matrix should be one less than V_0 along both dimentions. However, the
% user need not directly specify this, since any extra entries will just
% be dropped by the function anyway.
%
% FURTHER READING:
%
% Nagel, J. R, "Solving the Generalized Poisson Equation Using the
% Finite-Difference Method (FDM)" (Technical Report)
%
%==========================================================================
%
% James Nagel
% Department of Electrical and Computer Engineering
% University of Utah, Salt Lake City, Utah
% nageljr@ieee.org
% Copyright January 20, 2011
function V = Poisson_FDM_Solver_2D(V,BC,EPS,RHO,h)
% Permittivity of free-space (F/m).
EPS_0 = 8.854e-12;
% Extract simulation domain size.
[Ny,Nx] = size(V);
% Total number of unknowns to solve for.
L = Nx*Ny;
% Instantiate the dielectric coefficient matrices.
a0 = zeros(Ny,Nx);
a1 = zeros(Ny,Nx);
a2 = zeros(Ny,Nx);
a3 = zeros(Ny,Nx);
a4 = zeros(Ny,Nx);
% Short-hand index notation.
X1 = 2:Nx-1;
Y1 = 2:Ny-1;
X2 = 1:Nx-2;
Y2 = 1:Ny-2;
% Compute the dielectric coefficients.
a0(Y1,X1) = -( EPS(Y1,X1) + EPS(Y2,X1) + EPS(Y1,X2) + EPS(Y2,X2) );
a1(Y1,X1) = (0.5)*( EPS(Y1,X1) + EPS(Y2,X1) );
a2(Y1,X1) = (0.5)*( EPS(Y1,X2) + EPS(Y1,X1) );
a3(Y1,X1) = (0.5)*( EPS(Y2,X2) + EPS(Y1,X2) );
a4(Y1,X1) = (0.5)*( EPS(Y2,X1) + EPS(Y2,X2) );
% Separate coefficients into real and imaginary components. Matlab runs a
% lot quicker if we can limit ourselves to purely real-valued variable
% manipulation until the last moment.
a0r = real(a0); % Real
a1r = real(a1); %
a2r = real(a2); %
a3r = real(a3); %
a4r = real(a4); %
a0i = imag(a0); % Imaginary
a1i = imag(a1); %
a2i = imag(a2); %
a3i = imag(a3); %
a4i = imag(a4); %
% Check for charge densities.
if nargin == 3
b = zeros(Ny*Nx,1);
else
% Insert default h-value if not specified.
if nargin == 4
h = 1;
end
% Normalized charge matrix. This term comes from integrating the
% charge density over a square region with size h.
b = -RHO*h^2/EPS_0;
b = b(:); % Vectorize.
end
% Set the four corners to Dirichlet boundaries to avoid confusion in the
% algorithm. These points do not really matter anyway.
BC(1,1) = 1;
BC(1,Nx) = 1;
BC(Ny,Nx) = 1;
BC(Ny,1) = 1;
% Define flags for Neumann boundaries. Each case needs to be handled in
% its own unique way, so it helps to give each one a unique marker.
TOP_FLAG = -1;
BOTTOM_FLAG = -2;
LEFT_FLAG = -3;
RIGHT_FLAG = -4;
% Find the Neumann BCs on the edges.
Neumann_TOP = find(BC(1,:) ~= 1);
Neumann_BOTTOM = find(BC(Ny,:) ~= 1);
Neumann_LEFT = find(BC(:,1) ~= 1);
Neumann_RIGHT = find(BC(:,Nx) ~= 1);
% Set flags for Neumann boundaries.
BC(1,Neumann_TOP) = TOP_FLAG;
BC(Ny,Neumann_BOTTOM) = BOTTOM_FLAG;
BC(Neumann_LEFT,1) = LEFT_FLAG;
BC(Neumann_RIGHT,Nx) = RIGHT_FLAG;
% Initialize indices and values to fill the system matrix.
NZmax = 5*L; % Maximum possible number of nonzero elements.
I = zeros(NZmax,1); % i-indices.
J = zeros(NZmax,1); % j-indices.
Sr = zeros(NZmax,1); % Real Values.
Si = zeros(NZmax,1); % Imag Values.
idx = 1; % Nonzero entry index.
% Begin iterating over the unknowns and prepare to fill the system
% matrix. Filling the A-matrix is MUCH faster if you know the indices and
% values "a priori" for ALL non-zero elements. So rather than directly
% fill in the A-matrix, this loop stores all the nonzero (i,j) indices
% along with their values. The A-matrix is then instantiated directly
% from this information.
% Note that by convention, Matlab scans matrices columnwise when only
% a single element is indexed. So rather than waste time vectorizing the
% BC-matrix or initial-valued V-matrix, just use a single index to load
% and store values and remember this convention.
disp('Filling System Matrix');
for n = 1:L
% Check boundary condition.
switch BC(n)
% Dirichlet boundary.
case 1,
% A(n,n) = 1.
I(idx) = n;
J(idx) = n;
Sr(idx) = 1;
idx = idx + 1;
% Specify right-hand side as the voltage at this point.
b(n) = V(n);
% Top Neumann boundary.
case TOP_FLAG,
% A(n,n) = 1.
I(idx) = n;
J(idx) = n;
Sr(idx) = 1;
idx = idx + 1;
% A(n,n+1) = -1.
I(idx) = n;
J(idx) = n + 1;
Sr(idx) = -1;
idx = idx + 1;
% Bottom Neumann boundary.
case BOTTOM_FLAG,
% A(n,n) = 1.
I(idx) = n;
J(idx) = n;
Sr(idx) = 1;
idx = idx + 1;
% A(n,n-1) = -1.
I(idx) = n;
J(idx) = n - 1;
Sr(idx) = -1;
idx = idx + 1;
% Left Neumann boundary.
case LEFT_FLAG,
% A(n,n) = 1.
I(idx) = n;
J(idx) = n;
Sr(idx) = 1;
idx = idx + 1;
% A(n,n+Ny) = -1.
I(idx) = n;
J(idx) = n + Ny;
Sr(idx) = -1;
idx = idx + 1;
% Right Neumann boundary.
case RIGHT_FLAG,
% Set A(n,n) = 1
I(idx) = n;
J(idx) = n;
Sr(idx) = 1;
idx = idx + 1;
% A(n,n-Ny) = -1.
I(idx) = n;
J(idx) = n - Ny;
Sr(idx) = -1;
idx = idx + 1;
% Regular point. Apply the 5-point star.
otherwise,
% Convention for 5-point star:
%
% V2 | Indexing direction
% V3 V0 V1 |
% V4 \-/
%
% REMINDER: Single-valued indexing of a Matrix will scan
% COLUMN-WISE!
% V0 (center) term: A(n,n) = a0(n).
I(idx) = n;
J(idx) = n;
Sr(idx) = a0r(n); % Real
Si(idx) = a0i(n); % imaginary
idx = idx + 1;
% V1 (right) term: A(n,n+Ny) = a1(n)
I(idx) = n;
J(idx) = n + Ny;
Sr(idx) = a1r(n); % Real
Si(idx) = a1i(n); % imaginary
idx = idx + 1;
% V2 (top) term: A(n,n+1) = a2(n)
I(idx) = n;
J(idx) = n + 1;
Sr(idx) = a2r(n); % Real
Si(idx) = a2i(n); % imaginary
idx = idx + 1;
% V3 (left) term: A(n,n-Ny) = a3(n)
I(idx) = n;
J(idx) = n - Ny;
Sr(idx) = a3r(n); % Real
Si(idx) = a3i(n); % imaginary
idx = idx + 1;
% V4 (bottom) term: A(n,n-1) = a4(n)
I(idx) = n;
J(idx) = n - 1;
Sr(idx) = a4r(n); % Real
Si(idx) = a4i(n); % imaginary
idx = idx + 1;
end
end
% Throw out the leftover zeros.
I = I(1:idx-1);
J = J(1:idx-1);
Sr = Sr(1:idx-1);
Si = Si(1:idx-1);
% Combine real and imaginary coefficients into one vector.
S = Sr + 1j*Si;
% Fill the system matrix.
A = sparse(I,J,S,L,L,NZmax);
% Clear out memory before performing inversion. The matrix inversion
% process is HUGELY memory intensive, and will greedily hog up all the RAM
% it can get. Clearing out the major variables from the workspace will at
% least give a few extra MB of memory to this process. It won't be much,
% but every little bit helps if we're inverting a very large system.
clear I J S BC a0 a1 a2 a3 a4 a0r a0i Sr Si ;
clear a1r a1i a2r a2i a3r a3i a4r a4i Neumann_TOP Neumann_BOTTOM;
clear Neumann_LEFT Neumann_RIGHT EPS RHO;
% Invert the matrix. The slash operator is fantastic at doing this
% quickly for sparse matrices.
disp('Solving System Matrix');
V = A\b;
disp('Done!');
% Put the voltages back into a matrix.
V = reshape(V,Ny,Nx);
return;