Abnormal memory usage when working with large object collections

[foodaggression]

I'm calculating the magnetic field of a permanent magnet with a complex shape. To do this, I take the CAD drawing of the permanent magnet, export its surface as a STL Mesh and then use python to assemble a collection of magpylib triangle objects.

This process is fast, it only takes seconds to build a collection of 20k triangles - which is the lower bound for the meshing tolerance on the dimensions I care about. If I use a lower number of triangles, the field solution still changes significantly.

Calculating the field of this complex permanent magnet is also still relatively fast (around 10 seconds for 1000 points - testament to the excellence of magpylib!), but I run into problems when I want to calculate magnetic field maps at higher resolutions. Calling getB on this collection and a (100,100) array of locations uses more than 30GB of memory!

So what I currently do is building the magnet once, and then calling getB from within a loop along one axis of the map, and storing the result of each iteration by hand. This way, the field map completes in less than 2 minutes.

Is there a smarter way to do this? Does anybody know, why memory requirements of complex magnets explode with the number of observer points (I understand why the 20k triangles need quite a bit of memory)?

Any suggestions are appreciated!

[Alexboiboi]

Hi @foodaggression, nice to see that you make use of our latest Triangle class to build complex forms. At first glance it looks like you are pushing the limits of what magpylib is intended to be capable of, but it does not mean that there is no room for optimization ;)

Could you provide a self-contained example which shows the RAM issue, otherwise we would need to write one on our own to do some profiling?

Also, we are currently working on a new TriangularMesh class specifically designed to build complex body more easily and robustly. I don't know yet from the top of my head if it would help in your case but it may be worth have it a go.
Current work is progressing in the trimesh branch and some examples can be checked here. Note that the api is subjected to change, as this the feature is still under development.

Hope it helps in the meantime. We'll get back to you ;)

Sincerely, Alex -

[foodaggression]

I got around to testing the new TriangularMesh from the trimesh branch, and the problem seems to persist.

stl_mesh` = mesh.Mesh.from_file('example.stl')
magnet = magpy.magnet.TriangularMesh.from_triangular_facets(np.array([0,100,0]), facets=stl_mesh.vectors, validate_connected=False)

xs = np.linspace(-1e3,1e3,100)
ys = np.linspace(0,1e3,100)
grid = np.array([[(x,y,0) for x in xs] for y in ys])
B = magpy.getB(magnet, grid)

Nice work on the from_triangular_facets function, btw. It creates the 17k triangle magnet practically instantly! And it plays beautifully with data loaded by numpy-stl.

[Alexboiboi]

Thanks for your feedback ;)

When calculating the field, all the variables are tiled with the length of the flattened observer array and the number of triangles before feeding it to the core field function. This allows to vectorize the calculation and is almost always much faster than looping.
In the example you sent, the mesh object has 17k triangles and you want to compute 100*100 positions. The vertices parameters already gets very big in this case : 17000 * 100 * 100 * 3 * 8 bytes ≈ 3.8 GB. Considering that interim results with about the same size also need to be stored temporarily during computations it is not that surprising that 30GB+ (i got through with around 50GB) are needed.
Like you already pointed out, the alternative is to split-up the computation in smaller tasks and loop over them. In the end you can still apply a numpy summation over the relevant axis. It is up to you to choose to split the observers or the sources.

NB: There is a plan to switch the core computation to a lower, more efficient, lower level language. But this is not going to happen any time soon. You can track it here:

• MagLib - A Fortran-based core #579

[Alexboiboi]

@OrtnerMichael Feel free to chime in, if you have any better suggestion or any remarks ;)

[foodaggression]

Thank you so much for taking the time to explain all this!

One final bit of feedback: the new TriangularMesh consistently seems to give me ~20% slower performance when calling getB for the example object and 1000 points (when compared to just using a collection of 17k triangles forming the same object - as I did before). Results, of course, are identical.

I can submit example code if you're interested. Otherwise I consider this question answered.

[Alexboiboi]

Thank you so much for taking the time to explain all this!

No problem, glad to help ;)

One final bit of feedback: the new TriangularMesh consistently seems to give me ~20% slower performance when calling getB for the example object and 1000 points (when compared to just using a collection of 17k triangles forming the same object - as I did before). Results, of course, are identical.

This is expected since for the TriangularMesh source there is an additional check if the observers are inside or outside the body. The computation will be even slower if the observers are closer than the source bounding box, as the inside/outside check is much more computationally intensive (ray tracing). There is an option to opt out, but not available outside the core function, yet.

I can submit example code if you're interested. Otherwise I consider this question answered.

Sure ;)

[Alexboiboi]

Hi @foodaggression,

A way to get around the memory issue, is, as you mentioned, to divide the computation into chunks.

import magpylib as magpy
import matplotlib.pyplot as plt
import numpy as np
import pyvista as pv


def get_field_in_chunks(sources, points, nchunks=1, field="B", **kwargs):
    """Compute B-or-H-field for given sources and points in chunks. This allows to limit the memory
    usage of the fully vectorized getB or getH functions of the Magpylib library by dividing the
    observer points into subsets (`nchunks`). This should only be used if the number of sources AND
    the number of observers is very large, otherwise the computation time may decrease 
    significantly (e.g. TriangularMesh with 10000 facets + 1000 observers).
    
    Parameters
    ----------
    sources: source and collection objects or 1D list thereof
        Sources that generate the magnetic field. Can be a single source (or collection)
        or a 1D list of l source and/or collection objects.

    points: array_like
        Aarray_like positions of shape (n1, n2, ..., 3) where the field
        should be evaluated. All positions are given in units of [mm].
        
    nchunks: int
        Positive integer corresponding to the number of position chunks which the computation
        should be divided into.
        
    Notes
    -----
    See `magpylib.getB` or `magpylib.getH` for other keyword arguments.
    """
    
    Np = np.prod(points.shape[:-1]) # number of points
    nchunks = min(max(0, nchunks), Np) # make sure 0<nchunks<Np
    points = np.array(points, dtype=float)
    slices = [slice(i * Np // nchunks, (i + 1) * Np // nchunks) for i in range(nchunks)]
    func = getattr(magpy, f"get{field}")
    B = np.concatenate(
        [func(magnet, points.reshape(-1,3)[sl], **kwargs) for sl in slices]
    )
    return B.reshape(points.shape)


Nt = 10000 # target number of triangles
r = int(np.sqrt(Nt/2))
sphere  = pv.Sphere(theta_resolution=r, phi_resolution=r)
magnet = magpy.magnet.TriangularMesh.from_pyvista(
    magnetization=[0, 100, 0],
    polydata = sphere,
)

N = 100 # number points per side length
xs = np.linspace(-1e3, 1e3, N)
ys = np.linspace(0, 1e3, N)
grid = np.array([[(x, y, 0) for x in xs] for y in ys])

B = get_field_in_chunks(magnet, grid, nchunks=100)

This takes around 1min on my laptop.

We may integrate this into the main getB and getH function. Feel free to open a feature request if this is something that you would like to have.

Hope this helps ;)

Sincerely, Alex