# Maxwell’s Equations

<long unrelated intro>

I want to write a general Maxwell solver.

Energy and momentum are carried by electromagnetic fields through the Poynting vector and the energy density:

$$\vec{p}_{em}=\vec{S}/c^{2}=(1/\mu_{0})\vec{E}\times\vec{B}$$

$$u_{em}= (\epsilon_{0}/2)E^{2}+(1/2\mu_{0})B^{2}$$

They follow the following conservation laws:

$$\frac{\partial u_{em}}{\partial t}+\nabla\cdot\vec{S}+\vec{J}\cdot\vec{E}=0$$

$$\frac{\partial \vec{p}_{em}}{\partial t}-\nabla\cdot\sigma+\rho\vec{E}+\vec{J}\times\vec{B}=0$$

I was thinking it would be neat to have a finite volume method solver for Maxwell’s equations in terms of these conservation laws, but this is a flawed conception as the energy equations aren’t closed and IMO any way to close them ruins the neatness. So instead I decided to take the exact opposite extreme and just play around with a finite difference scheme.

The electric and magnetic fields can be written in terms of the potentials $$\phi$$ and $$\vec{A}$$. In the Lorenz gauge, Maxwell’s equations are:

$$\nabla^{2}\phi-1/c^{2}\frac{\partial^{2} \phi}{\partial t^{2}}=-\rho/\epsilon_{0}$$

$$\nabla^{2}\vec{A}-1/c^{2}\frac{\partial^{2} \vec{A}}{\partial t^{2}}=-\mu_{0}\vec{J}$$

Which is the inhomogenous wave equation four times. I found a stencil for the 3D Laplacian in [R.C. O’Reilly and J.M. Beck Int. J. Numer. Meth. Engng 2006; 00-1-16] which also conducts stability analysis for solving the wave equation and heat equations using a forward Euler scheme.

$$\nabla^{2}\phi\approx\frac{3}{13h^{2}}\Big(\sum_{j\in\mathcal{N_{f}}} \phi_{j}+\frac{1}{2} \sum_{j\in\mathcal{N}_{e}} \phi_{j}+\frac{1}{3} \sum_{j\in\mathcal{N}_{c}} \phi_{j} – \frac{44}{3}\phi_{i}\Big)$$

Just as an example in C++:


#define n_points 512
#define tot_points n_points*n_points*n_points
#define dx 0.1
#define dt 0.01
#define cec -(44.0/3.0)*(3.0/(13.0*dx*dx)) //center,face,..
#define fac (3.0/(13.0*dx*dx))
#define edc 0.5*(3.0/(13.0*dx*dx))
#define coc (1.0/3.0)*(3.0/(13.0*dx*dx))
#define forxyz(X) for(int x = 1; x &lt; n_points-1; x++) {for(int y = 1; y &lt; n_points-1; y++) {for(int z = 1; z &lt; n_points-1; z++) { X }}}

int c(int x, int y, int z) {
return n_points*n_points*x+n_points*y+z;
}

void laplacian(std::array&lt;float,tot_points> *arrout, std::array&lt;float,tot_points> *arrin) {
(*arrout) = {};//1 sec out of 15
//9 3 slices

float temp;

for(int i = -1; i &lt; 2; i++) {
for(int j = -1; j &lt; 2; j++) {
int c1;
int c2;

if((i==0)and(j==0)) {//face center face
forxyz(
c1=c(x+i,y+j,z);
c2=c(x,y,z);
(*arrout)[c2]+=fac*(*arrin)[c1-1]+cec*(*arrin)[c1]+fac*(*arrin)[c1+1];
);

} else if((abs(i)==1)and(abs(j)==1)) { //corner edge corner
forxyz(
c1=c(x+i,y+j,z);
c2=c(x,y,z);
(*arrout)[c2]+=coc*(*arrin)[c1-1]+edc*(*arrin)[c1]+coc*(*arrin)[c1+1];
);

} else {//edge face edge
forxyz(
c1=c(x+i,y+j,z);
c2=c(x,y,z);
(*arrout)[c2]+=edc*(*arrin)[c1-1]+fac*(*arrin)[c1]+edc*(*arrin)[c1+1];
);
}
}
}
return;
}



To be updated with CUDA, method of relaxation, …

~-~-~-~-~-~-~-~