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:


\(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                                                                                                                                                

       } else if((abs(i)==1)and(abs(j)==1)) { //corner edge corner                                                                                                                            

       } else {//edge face edge                                                                                                                                                                

To be updated with CUDA, method of relaxation, …

Author: Garth Whelan


Leave a Reply

Your email address will not be published. Required fields are marked *