<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.

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 < n_points-1; x++) {for(int y = 1; y < n_points-1; y++) {for(int z = 1; z < 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<float,tot_points> *arrout, std::array<float,tot_points> *arrin) {

(*arrout) = {};//1 sec out of 15

//9 3 slices

float temp;

for(int i = -1; i < 2; i++) {

for(int j = -1; j < 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, …