Solvers and Performance#
apply_demag offers two solvers that share the same interaction model (the
analytical volume-averaged Newell tensor for parallel cuboid cells, point
matching otherwise) and therefore agree to solver tolerance for any input:
solver="direct"(default) — assembles the dense system matrix and solves it exactly with LAPACK. Memory grows as (3N)² and time as N³, so it is the right choice up to a few thousand cells.solver="iterative"— solves the same system matrix-free with GMRES. Cells on uniform grids (the output ofmesh_Cuboid) are handled by an O(N log N) FFT convolution, so meshes with tens of thousands of cells stay fast and memory-light. If GMRES cannot reachsolver_tolwithinmax_iteriterations, aRuntimeErroris raised — a partially converged result is never returned silently.
This example verifies the agreement on typical use cases and measures the performance crossover.
import time
import tracemalloc
import magpylib as magpy
import numpy as np
import pandas as pd
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from magpylib_material_response.demag import apply_demag
from magpylib_material_response.meshing import mesh_Cuboid
magpy.defaults.display.backend = "plotly"
def max_rel_diff(coll_a, coll_b):
"""Maximum relative polarization difference between two solved collections."""
pol_a = np.array([s.polarization for s in coll_a.sources_all])
pol_b = np.array([s.polarization for s in coll_b.sources_all])
return np.abs(pol_a - pol_b).max() / np.abs(pol_a).max()
Both solvers agree — typical use cases#
Three configurations that exercise all interaction paths: a single meshed magnet (FFT grid), several bodies with different cell sizes plus a standalone magnet (grid + cross + generic blocks), and a rotated mesh with anisotropic susceptibility (rotation handling).
# 1. single meshed cuboid magnet (SI units: m, T)
cube = magpy.magnet.Cuboid(polarization=(0, 0, 1), dimension=(1e-3, 1e-3, 1e-3))
case_single = magpy.Collection(mesh_Cuboid(cube, target_elems=343))
# 2. two meshed bodies with different cell sizes + one standalone magnet
body_a = magpy.magnet.Cuboid(polarization=(0, 0, 1), dimension=(1e-3, 1e-3, 1e-3))
body_b = magpy.magnet.Cuboid(
polarization=(1, 0, 0), dimension=(2e-3, 1e-3, 1e-3), position=(2e-3, 0, 0)
)
lone = magpy.magnet.Cuboid(
polarization=(0, 0, 1), dimension=(0.5e-3, 0.5e-3, 0.5e-3), position=(0, 2e-3, 0)
)
case_mixed = magpy.Collection(mesh_Cuboid(body_a, 125), mesh_Cuboid(body_b, 125), lone)
# 3. rotated meshed cuboid with anisotropic susceptibility
cube_rot = magpy.magnet.Cuboid(polarization=(0, 0, 1), dimension=(1e-3, 1e-3, 1e-3))
cube_rot.rotate_from_angax(35, (1, 1, 0))
case_rotated = magpy.Collection(mesh_Cuboid(cube_rot, 125))
magpy.show(
{"objects": case_single, "col": 1},
{"objects": case_mixed, "col": 2},
{"objects": case_rotated, "col": 3},
)
cases = {
"single meshed magnet": (case_single, 3.0),
"two bodies + standalone": (case_mixed, 3.0),
"rotated, anisotropic": (case_rotated, (0.3, 0.1, 0.5)),
}
rows = []
for label, (coll, sus) in cases.items():
n = len(coll.sources_all)
sus_list = [sus] * n
coll_direct = apply_demag(coll, susceptibility=sus_list, solver="direct")
coll_iter = apply_demag(
coll, susceptibility=sus_list, solver="iterative", solver_tol=1e-8
)
rows.append(
{
"use case": label,
"cells": n,
"max rel. difference": f"{max_rel_diff(coll_direct, coll_iter):.1e}",
}
)
pd.DataFrame(rows)
| use case | cells | max rel. difference | |
|---|---|---|---|
| 0 | single meshed magnet | 343 | 1.6e-08 |
| 1 | two bodies + standalone | 254 | 7.6e-09 |
| 2 | rotated, anisotropic | 125 | 1.4e-09 |
The two solvers match to the requested solver_tol in every configuration — the
choice between them is purely a performance trade-off.
Performance scaling#
Wall time and peak memory for a single meshed cuboid magnet of growing cell
count. Absolute seconds depend heavily on the machine, so times are reported
relative to the direct solve of the smallest mesh — the ratios are much more
portable. Peak memory is measured with tracemalloc (which tracks NumPy
allocations) in a separate run, since tracing skews wall time.
def measure_solve(target_elems, solver):
cube = magpy.magnet.Cuboid(polarization=(0, 0, 1), dimension=(1e-3, 1e-3, 1e-3))
coll = magpy.Collection(mesh_Cuboid(cube, target_elems=target_elems))
n = len(coll.sources_all)
t0 = time.perf_counter()
apply_demag(coll, susceptibility=3.0, solver=solver, solver_tol=1e-8)
dt = time.perf_counter() - t0
# separate traced run (tracing skews wall time); apply_demag works on a
# copy, so the same collection can be reused
tracemalloc.start()
apply_demag(coll, susceptibility=3.0, solver=solver, solver_tol=1e-8)
_, peak = tracemalloc.get_traced_memory()
tracemalloc.stop()
return {"solver": solver, "cells": n, "time": dt, "peak_mb": peak / 1e6}
# warm-up both solvers once, so one-time costs (imports, BLAS/FFT
# initialisation) do not leak into the first timed measurement
warmup = magpy.Collection(
mesh_Cuboid(magpy.magnet.Cuboid(polarization=(0, 0, 1), dimension=(1e-3, 1e-3, 1e-3)), 64)
)
apply_demag(warmup, susceptibility=3.0, solver="direct")
apply_demag(warmup, susceptibility=3.0, solver="iterative", solver_tol=1e-8)
sizes_direct = [216, 512, 1000, 1728]
sizes_iterative = [216, 512, 1000, 1728, 4096, 8000]
records = [measure_solve(s, "direct") for s in sizes_direct]
records += [measure_solve(s, "iterative") for s in sizes_iterative]
perf_df = pd.DataFrame(records)
n_ref = int(perf_df["cells"].min())
t_ref = perf_df[(perf_df["solver"] == "direct") & (perf_df["cells"] == n_ref)][
"time"
].iloc[0]
perf_df["rel_time"] = perf_df["time"] / t_ref
perf_df.pivot(index="cells", columns="solver", values=["rel_time", "peak_mb"]).round(2)
| rel_time | peak_mb | |||
|---|---|---|---|---|
| solver | direct | iterative | direct | iterative |
| cells | ||||
| 216 | 1.00 | 0.68 | 14.87 | 0.79 |
| 512 | 4.67 | 2.26 | 82.51 | 1.86 |
| 1000 | 18.17 | 2.28 | 313.39 | 3.48 |
| 1728 | 60.97 | 3.76 | 934.04 | 6.01 |
| 4096 | NaN | 9.53 | NaN | 14.17 |
| 8000 | NaN | 18.52 | NaN | 27.56 |
series_style = {
"direct": {"color": "#2a78d6", "symbol": "circle"},
"iterative": {"color": "#1baf7a", "symbol": "square"},
}
fig = make_subplots(
rows=1,
cols=2,
subplot_titles=("wall time (relative)", "peak memory"),
horizontal_spacing=0.12,
)
for solver, style in series_style.items():
df = perf_df[perf_df["solver"] == solver]
line = {"color": style["color"], "width": 2}
marker = {"color": style["color"], "symbol": style["symbol"], "size": 9}
fig.add_trace(
go.Scatter(
x=df["cells"],
y=df["rel_time"],
name=solver,
mode="lines+markers",
line=line,
marker=marker,
hovertemplate="%{x} cells: %{y:.1f}×<extra>" + solver + "</extra>",
),
row=1,
col=1,
)
fig.add_trace(
go.Scatter(
x=df["cells"],
y=df["peak_mb"],
name=solver,
showlegend=False,
mode="lines+markers",
line=line,
marker=marker,
hovertemplate="%{x} cells: %{y:.0f} MB<extra>" + solver + "</extra>",
),
row=1,
col=2,
)
fig.add_annotation(
x=np.log10(df["cells"].iloc[-1]),
y=np.log10(df["rel_time"].iloc[-1]),
text=solver,
font={"color": style["color"]},
xanchor="left",
xshift=12,
showarrow=False,
xref="x",
yref="y",
)
fig.update_xaxes(title_text="number of cells", type="log")
fig.update_yaxes(type="log")
fig.update_yaxes(title_text=f"wall time (× direct at {n_ref} cells)", row=1, col=1)
fig.update_yaxes(title_text="peak memory (MB)", row=1, col=2)
fig.update_layout(
title="apply_demag scaling — single meshed cuboid",
legend={"orientation": "h", "yanchor": "bottom", "y": 1.08},
template="plotly_white",
)
fig.show()
The direct solver’s N³ time slope and (3N)² memory slope take over in the low thousands of cells, while the FFT-accelerated iterative solver stays almost flat in both panels. Beyond the crossover the gap widens rapidly — the dense matrix becomes the hard limit (extrapolating the right panel, N = 27 000 would already need a ~50 GB matrix, while the iterative solver handles it in seconds within a few hundred MB).
Model topology matters#
The solvers assemble the interaction operator from structure clusters (uniform
grids of identical parallel cells → FFT / analytical blocks, everything else →
point-matched magpy.getH), so the model topology decides which paths do the
work — and how much the solver choice matters. Here the same comparison runs on
characteristic topologies of similar total cell count (single-run timings —
indicative, not statistics).
from collections import Counter
from magpylib_material_response.demag_fft import analyze_collection
from magpylib_material_response.meshing import mesh_Cylinder
def structure_kinds(coll):
"""Summarize the structure clusters the solvers will work with."""
_, clusters = analyze_collection(coll.sources_all)
kinds = Counter(c["kind"] for c in clusters)
return " + ".join(f"{v} {k}" for k, v in sorted(kinds.items()))
def bench_topology(label, coll):
n = len(coll.sources_all)
sus = [0.5] * n
t0 = time.perf_counter()
coll_direct = apply_demag(coll, susceptibility=sus, solver="direct")
t_direct = time.perf_counter() - t0
t0 = time.perf_counter()
coll_iter = apply_demag(coll, susceptibility=sus, solver="iterative", solver_tol=1e-8)
t_iter = time.perf_counter() - t0
return {
"topology": label,
"cells": n,
"structure": structure_kinds(coll),
"iter / direct time": round(t_iter / t_direct, 2),
"agreement": f"{max_rel_diff(coll_direct, coll_iter):.0e}",
}
def block(pos=(0, 0, 0)):
return magpy.magnet.Cuboid(
polarization=(0, 0, 1), dimension=(1e-3, 1e-3, 1e-3), position=pos
)
def grid_positions(k, pitch=1.8e-3):
side = int(np.ceil(np.sqrt(k)))
return [(i * pitch, j * pitch, 0) for i in range(side) for j in range(side)][:k]
topologies = {
"single body": magpy.Collection(mesh_Cuboid(block(), 1000)),
"4 bodies": magpy.Collection(
*(mesh_Cuboid(block(p), 250) for p in grid_positions(4))
),
"24 bodies": magpy.Collection(
*(mesh_Cuboid(block(p), 48) for p in grid_positions(24))
),
}
meshes = []
for i, p in enumerate(grid_positions(24, pitch=2.2e-3)):
m = mesh_Cuboid(block(p), 48)
m.rotate_from_angax(i * 15.0, (0, 0, 1), anchor=p)
meshes.append(m)
topologies["24 bodies, individually rotated"] = magpy.Collection(*meshes)
cyl = magpy.magnet.Cylinder(polarization=(0, 0, 1), dimension=(1e-3, 1e-3))
topologies["meshed cylinder"] = magpy.Collection(mesh_Cylinder(cyl, 300))
for label, coll in topologies.items():
coll.style.label = label
magpy.show(
*(
{"objects": coll, "row": i // 3 + 1, "col": i % 3 + 1}
for i, coll in enumerate(topologies.values())
),
)
topo_rows = [
bench_topology(label, coll)
for label, coll in topologies.items()
if label != "meshed cylinder"
]
topo_rows.append(bench_topology("meshed cylinder", topologies["meshed cylinder"]))
pd.DataFrame(topo_rows)
| topology | cells | structure | iter / direct time | agreement | |
|---|---|---|---|---|---|
| 0 | single body | 1000 | 1 grid | 0.16 | 4e-09 |
| 1 | 4 bodies | 1008 | 4 grid | 0.44 | 6e-09 |
| 2 | 24 bodies | 1152 | 24 grid | 0.52 | 2e-08 |
| 3 | 24 bodies, individually rotated | 1152 | 24 grid | 0.93 | 2e-08 |
| 4 | meshed cylinder | 350 | 1 generic | 1.00 | 3e-09 |
What the rows show:
Single body / a few parallel bodies — everything runs on the FFT and analytical Newell paths; the iterative solver is ahead and scales far better (previous section).
Dozens of bodies — the pairwise cross-blocks (one per body pair) start to dominate the operator build. The iterative solver still wins, but its advantage shrinks as the body count grows at fixed total cells.
Bodies rotated individually — cross-blocks between differently oriented bodies fall back to point matching, which both solvers share: expect near parity. The self-blocks of each body still use their own rotated FFT grid.
Meshed cylinder — non-cuboid cells (here Cylinder and CylinderSegment, likewise tetrahedra from
mesh_TriangularMesh— see the tetrahedral mesh example) take the point-matched generic path. Both solvers spend nearly all time evaluating the cells’ analytical fields, so the solver choice barely matters — the cell type is the cost driver. Prefer cuboid meshes when the geometry allows it.
Choosing a solver#
Model topology |
Interaction paths |
Recommendation |
|---|---|---|
One meshed body, up to a few thousand cells (~2000–3000) |
FFT / analytical |
either; |
One meshed body, large |
FFT self-block |
|
Few parallel bodies |
FFT + analytical cross-blocks |
|
Dozens of bodies |
cross-blocks dominate the build |
|
Bodies rotated differently |
point-matched cross-blocks |
either — build cost is shared, expect parity |
Non-cuboid cells (cylinder, tetrahedra) |
point-matched, dense |
either — cell field evaluation dominates; keep counts moderate |
Two knobs control the iterative solver: solver_tol (relative residual, default
1e-6) and max_iter (default 50). Non-convergence raises a RuntimeError
with guidance rather than returning an inaccurate result. To see the detected
structure for your own model, enable the package logging
(magpylib_material_response.configure_logging()) — the cluster summary is
logged at the start of every solve.