leon-vv / traceon Goto Github PK

In commits fa32fe1 and 214cef2 we have lost some functionality to be able to push out v0.1 sooner. We should consider bringing back some of this functionality (especially automatic focus finding and fitting of aberration coefficients) behind a new shiny API.

We know that the Maxwell equations are linear. This means that if we solve for every electrode (unity excitation) we can use superposition to find the fields at different excitation. In terms of computation this means we run the solver N terms, where N is the number of electrodes.

Now, when we use a direct solver, according to the (Numpy documentation)[https://numpy.org/doc/stable/reference/generated/numpy.linalg.solve.html] behind the scenes uses a LAPACK routine which computes the LU decomposition. It would therefore be faster to store this LU decomposition after solving for the first electrode and reuse this decomposition for the other electrodes. In other words, we solve the system with multiple right hand sides, which is computationally cheaper.

This might have a big impact on the exposed API. We have to think through how to expose the superposition of fields in the API in a neat way. We need to be able to rescale the fields after a solution step.

Ideas to improve the method:

• $A^{\nu+1}_m$ should be the second derivative (in z direction) of coefficient $A^{\nu}_m$. We could use this to improve the accuracy of the cubic spline interpolation along the z axis.
• The current theta interpolation uses values computed using numerical integration in Python using 64 bit floating point numbers. Therefore they can be expected to be accurate to at most a small multiple of the 64 bit float machine epsilon. Since for the pre-computation we can expend an arbitrary amount of computing power, we need to make sure these values are correct down to the max precision achievable with 64 bit floating point numbers. We can achieve this using arbitrary precision floating point number numerical integration.
• The current theta interpolation uses a cubic spline interpolation. We can consider other interpolation types, also we can make the cubic spline interpolation more accurate by explicitly providing the derivatives, which we can compute to arbitrary precision ourselves.

Dielectric materials can readily be added to the simulations if we allow for Dirichlet boundary conditions. In other words conditions on the electric field instead of the potential. If I understand correctly, the electric field at the boundary of a dielectric and the bound surface charge induced there are related by

we therefore need to be able to generate normal vectors at the boundaries. This should be possible with GMSH getNormal function.

We use STEP_MIN and STEP_MAX to guide the particle tracing. However, due to the adaptive step size control do we really need these values? Some more research is needed.

The special geometries included in geometry.py should be in their own benchmark files. Not in the geometry module.

Could be hosted in the meantime on my personal webpage.

Some tests/examples on how to use Geometry.MEMSStack

The interpolation technique used to speed up ray tracing can be improved

It seems like matplotlib 3D renderer is a bit unreliable. Maybe switching to something like Mayavi with a proper 3D renderer would be beneficial.

We have to get rid of the Nmax variable in tracing.py. The particle should be traced indefinitely.

We should clearly specify which symmetries are supported and for which excitation (electrostatic/magnetostatic).

For finite length with planar symmetry we have to neglect edge effects as the problem would otherwise become solvable only in 3D. One could image not implementing the infinite length case as it can be approximated easily by increasing the finite length to a large value (I believe this is the approach Comsol takes).

Full support for a symmetry implies the following:

• Ability to solve for relevant fields (and compute axial derivatives)
• Ability to employ smart (Taylor) interpolation to do very fast electron tracing
• Comparison with either literature values or competing software suites showing high accuracy
• Documentation of the relevant functions and formulas

When simulating an entire microscope, we need to be able to trace an electron through multiple elements. The most obvious approach to achieve this would be to define a new SystemTracer which accepts multiple fields. Translation and rotations of the fields inside the 'global' coordinate system should be supported.

We need to update the scripts to the new API

Currently the only tests of the software suite are the validations in /validation.

We need a bigger test suite that clearly distinguishes between the following classes of tests:

• Accuracy test, comparisons with either literature values or alternative software showing the accuracy of the solvers
• Performance test, clearly identifying the accuracy/time tradeoff of the solvers
• Later: API tests, showing the library updates do not break backwards compatibility (only useful after v1.0)

Each class of test should use the same function to display the results.

Currently we use Simpson integration rule (see solver/util.py). It might be beneficent to instead use a Gaussian quadrature. We need to explore the performance/accuracy benefits. Hopefully using the better tests of issue #4.

Currently the numerical charge values are expressed as $\sigma / \pi \epsilon_0$. I think it would be semantically clearer if we expressed the quantity as $\sigma / \epsilon_0$. The reason is that the $\pi$ value only makes sense in cylindrical symmetry, but we need to implement more symmetries.

I remember reading somewhere (forgot the source) that the accuracy of the BEM is optimal when the charges on the elements are roughly equal. This means that once we find the solution for a given mesh size, we can refine the mesh in an intelligent way by considering which line elements have the most charge.

It looks like using multiple threads doesn't always speed up the program (observed especially in 3D). We have to investigate what causes these performances issues

For axial symmetry we have an elegant interpolation technique relying on high order axial derivatives of the potential. However, the computation time can quickly be dominated by the computation of the derivatives. Also, for the 3D and planar symmetric cases it's not obvious the computation of the high order derivatives will be equally elegant.

For this reason we need a new, reliable interpolation technique which is not as cumbersome to set up. Hermite interpolation seems the most obvious.

The Numba JIT compile times are getting out of hand. Also passing a function to a Numba compiled function (like we do for line_integral) is not properly supported in Numba with regards to caching. This means that those functions need to be recompiled often...

A proper C or Cython or Rust module would probably fix these issues.

For large 3D problems the direct solver can take a long time. Maybe good idea to support an iterative solver.

We currently use a axial Taylor expansion to quickly calculate the necessary fields. We need a script that for any geometry can plot the accuracy of the interpolation along the axis (for different distances from the axis).

Quoting from Peter W. Hawkes. Pinciples of Electron Optics, Basic Geometrical Optics. 2017.

Section 10.2.1:

Another advantage is the possibility of using more than the minimum number N of trial
functions and variables Q k . A favourable choice is, for instance, the use of function values
Y(s n ) and derivatives Y’(s n ) at the nodes, hence M = 2N. In each boundary element the
function Y(s) can then be approximated by a cubic Hermite polynomial, which is fairly
accurate. Each node has two degrees of freedom, which gives a better result than the simple
one-degree approximation.

After solving for the field, currently an ugly tuple is returned containing all the relevant data. It would be better to return a Field class containing a reference to the excitation and the geometry.

Now sometimes we over sample by a lot, for example validation/preikszas.py spends way too much time calculating the required derivatives. Maybe a good first guess would be to take the sampling equal to the mesh size.

Currently we only take the analytical derivative up to order 4 in solver/__init__.py function get_axial_derivatives. After using the assumption r_0 = 0 to simplify the formulas (see #3) we should be able to calculate the derivatives of higher order analytically.

See also formula 13.16a in Peter W. Hawkes. Principles of Electron Optics, Basic Geometrical Optics. 2017.

See here: https://about.zenodo.org/

We need to have a validation for this case. Convergence to the infinite case when the finite length grows large would be acceptable for a first validation.

The relevant formula for the potential seems to be given by this page.

Don't only return the position of the particles, but also the corresponding time steps

Many function docstrings are now obsolete

We currently use N as input to approximate the number of mesh elements. However, the approximation is really bad. Better would be to simply take the mesh size as input.

Currently in traceon/solver/radial_symmetry.py we have not used the assumption that r_0 = 0. Using this assumption greatly simplifies the formulas and therefore likely speeds up the calculation of the derivatives.

The particle mass and charge should be adjustable.

The validations act as examples too.. but are not very penetrable. We need better example files, for example a simple einzel lens and a 3D deflector ?

Currently we compute the self voltage (voltage produced at an element by the charge on the same element) by making sure the integration doesn't step into the singularity at the center. It would be better to figure out a more accurate procedure to calculate the self voltage. I recall that a exact formula was claimed for the triangle case in:

Computational Charged Particle Optics. Thomas Verduin. 2013.