Enhancement: Implementation of FastMultipole method also for magnetostatic problems.

[renereimann]

Currently magnetostatic problems are either solved using zonal harmonics (if the geometry is cylinder symmetric) or has to be solved using numerical integration of the Biot-Savart law (for all non cylinder symmetric setups). For sufficient complex current / wire distributions the numerical integration of the Biot-Savart law is a major bottleneck at run time as for each particle step the magnetic field has to be reevaluated including the numerical integration over all wires.

For electrostatic fields the FastMultipole method was implemented that allows to precompute boundary charges on a surface and using a tree representation for fast evaluation at run time. The method is based on the fact that the volume of interest is charge free and thus the E field is divergent and curl free and thus can be described by a potential that has to fulfill the Laplace equation.

For magnetic fields we usually are also interested in the current free region (we are not tracking particles through wires) and thus the magnetic field is divergent and curl free and can also be described by a potential that has to fulfill the Laplace equation.

Thus we can make use of the same mathematical description as in the case of electromagnetic problems which are described by a set of surface charges sigma_i on an enclosing surface. For electrostatic fields the surface charges have to be derived from Direchlet and Neuman boundary conditions. The surface charges are given by sigma(r') = - eps_0 normal(surface) * grad(potential(r')) .
We can make use of -grad(potential(r')) = B(r') and thus can directly calculate sigma(r') = const * normal(surface) * B(r'). In this case we can calculate B(r') using the Biot-Savart law at initialization time to get the surface charges and the FastMultipole representation. This has a speed up because 1. we have to evaluate B(r') only at a surface and not for the full volume, 2. we only do it at initialization time once and not at every particle step and 3. we can save the representation for other executions.

Implementation: We need a new class for magnetostatic problems what calculates "pseudo" surface charges using the Biot-Savart integration (which is already implemented). Based on the "pseudo" surface charges we generate the tree representation of the FastMultipole method (already implemented) and at runtime we use the FastMultipole method (already implemented) also for magnetic fields. Thus implementation is basically the work of writing wrapper functions to already existing tools.

[renereimann]

Currently the runtime for an Ioffe style magnetic trap takes me >10h for a single electron. I will work on the implementation over the next weeks but if others feel like they are faster please join the effort.

[zykure]

Hi René, that is a very interesting suggestion! And it sounds even better when you make clear that most of the necessary tools are already in place, just not for the magnetostatic fields.

At KATRIN usually the electric field is the bottleneck, so there has not been much pressure to optimize the magnetic fields as well. (And most magnets are symmetric solenoids / air coils anyway, so we can use the ZH approximation.) Hopefully with this new method you can cut down the simulation time to well below the 10h you currently need. Good luck!

[jpbarrett]

Hi Rene, Jan,

It's been a while since I've looked at the fast multipole code, but I'd be happy to answer questions about it if you go ahead with adding a magnetostatic solver. A long while ago I had considered implementing something like this, but I got busy with other projects and it never got very far. Originally, I was thinking of doing this with a cartesian 3-vector of multipole coefficients (along the lines of this paper: paper1). Each component of the vector would just transform independently as a separate scalar in the same way as the original (electrostatic) field map, so the only real differences would be: (1) managing the fast-multipole process 3 times, (2) calculating the vector multipole coefficients from the current distributions, and then (3) computing the field from the vector local coefficients.

I think the magnetic scalar potential approach you suggest could work too, you'd just need to be careful about distinguishing between regions with and without free currents. I think this approach is commonly done for magnetic material simulations where the volume elements are cubic (e.g. paper2), and there are no free currents involved, just polarized material. In that particular case the surfaces over which the magnetic 'charges' are glued are pretty straightforward since they are just the surface of each cubic element of the magnetic material.

In the case where you do have free currents on wires this might be more complicated, I suppose the work-around would be to define and discretize a surface enclosing your region of interest and then compute the appropriate equivalent magnetic charges over that surface from the global magnetic field and/or current sources. The exact details of how to do that discretization to get the required accuracy might be complicated though, but should be less of a problem if you can guarantee that the region of interest where your particles live is well separated from the 'magnetic charge surface'.