Microcubed vs. Ubermag/OOMMF: stray field and gradients¶
This notebook compares the analytical cuboid solution in microcubed with an Ubermag/OOMMF airbox calculation. It follows the setup of the Ubermag stray-field tutorial.
For an apples-to-apples comparison the magnetisation is kept uniform: no relaxation is performed. OOMMF returns \(H\) in A/m, which is converted to \(B=\mu_0H\) in tesla. Microcubed uses nanometres for geometry in this notebook, so its analytical gradient is in T/nm. OOMMF gradients use second-order finite differences on the same 5 nm grid (central in the interior, one-sided at the airbox edges).
import os
from pathlib import Path
from tempfile import TemporaryDirectory
import discretisedfield as df
import matplotlib.pyplot as plt
import micromagneticmodel as mm
import numpy as np
import oommfc as oc
import pandas as pd
from microcubed import Magnet
plt.style.use("default")
AXIS_NAMES = ("x", "y", "z")
Optional installation and OOMMF runner¶
Install with uv sync --extra examples --extra comparison. The comparison
extra is not required by Microcubed, the core examples, or ordinary docs builds.
OOMMF itself is a separate solver: install it using the
Ubermag instructions.
Set OOMMFTCL=/absolute/path/to/oommf.tcl for a Tcl installation; otherwise
use oommfc’s configured runner (for example its Docker runner).
Run uv run --extra examples --extra comparison python tools/notebooks.py --check --execute --include-optional.
To execute this page within Sphinx, also set MICROCUBED_RUN_OOMMF=1 and install the docs extra.
Ordinary docs builds display the source without running this optional solver.
# Runner discovery can perform a small validation calculation. Keep its files temporary.
with TemporaryDirectory(prefix="microcubed-oommf-runner-") as runner_directory:
previous_directory = Path.cwd()
try:
os.chdir(runner_directory)
if os.environ.get("OOMMFTCL"):
oc.runner.runner = oc.oommf.TclOOMMFRunner(os.environ["OOMMFTCL"])
runner = oc.runner.runner
finally:
os.chdir(previous_directory)
print(f"OOMMF runner ready: {runner}")
Ubermag/OOMMF stray field¶
Only the demagnetising effective field is computed. The temporary OOMMF files are removed after the returned field has been loaded.
with TemporaryDirectory(prefix="microcubed-oommf-") as output_directory:
H_ubermag = oc.compute(
system.energy.demag.effective_field,
system,
dirname=output_directory,
append=False,
verbose=1,
)
B_ubermag = MU0 * H_ubermag.array # shape: (nx, ny, nz, component), T
assert B_ubermag.shape == (*mesh.n, 3)
assert np.isfinite(B_ubermag).all()
Microcubed field at the identical cell centres¶
OOMMF coordinates are converted from metres to nanometres before being passed to microcubed. Values inside the magnet are intentionally masked by microcubed and excluded from the comparison.
x, y, z = (np.asarray(axis_cells) for axis_cells in mesh.cells)
xx, yy, zz = np.meshgrid(x, y, z, indexing="ij")
coordinates_m = np.stack([xx, yy, zz], axis=-1)
points_nm = np.moveaxis(coordinates_m, -1, 0).reshape(3, -1) * 1e9
cube = Magnet(
size=[100, 100, 100],
center=[0, 0, 0],
magnetization=[0, 0, MS],
)
B_microcubed = cube.Bfield(points_nm).T.reshape(*mesh.n, 3)
outside = np.any(np.abs(coordinates_m) > MAGNET_HALF_SIZE, axis=-1)
distance_to_cube = np.linalg.norm(np.maximum(np.abs(coordinates_m) - MAGNET_HALF_SIZE, 0), axis=-1)
# Keep two cells away from the discontinuous material boundary for global metrics.
field_mask = outside & (distance_to_cube >= 2 * CELL) & np.all(np.isfinite(B_microcubed), axis=-1)
print(f"Compared exterior field cells: {field_mask.sum():,}")
def error_metrics(reference, candidate, mask=None):
reference = np.asarray(reference)
candidate = np.asarray(candidate)
valid = np.isfinite(reference) & np.isfinite(candidate)
if mask is not None:
valid &= mask
reference = reference[valid]
candidate = candidate[valid]
difference = candidate - reference
rmse = np.sqrt(np.mean(difference**2))
reference_rms = np.sqrt(np.mean(reference**2))
return {
"RMSE": rmse,
"NRMSE": rmse / reference_rms if reference_rms else np.nan,
"max_abs_error": np.max(np.abs(difference)),
"correlation": np.corrcoef(reference, candidate)[0, 1],
}
field_metrics = pd.DataFrame(
{
f"B{axis}": error_metrics(B_ubermag[..., component], B_microcubed[..., component], field_mask)
for component, axis in enumerate(AXIS_NAMES)
}
).T
assert (field_metrics["NRMSE"] < 0.01).all(), field_metrics
field_metrics
Field maps on an exterior plane¶
Each column is one field component. Rows show microcubed, Ubermag/OOMMF, and their signed difference.
plane_index = int(np.argmin(np.abs(z + 72.5e-9)))
plane_z_nm = z[plane_index] * 1e9
x_nm, y_nm = x * 1e9, y * 1e9
def comparison_figure(reference, candidate, quantity_names, title, unit):
difference = candidate - reference
fig, axes = plt.subplots(3, len(quantity_names), figsize=(15, 11), constrained_layout=True)
rows = ((candidate, "microcubed"), (reference, "Ubermag/OOMMF"), (difference, "difference"))
for column, quantity in enumerate(quantity_names):
for row, (data, row_name) in enumerate(rows):
values = data[..., column]
limit = (
max(np.nanmax(np.abs(reference[..., column])), np.nanmax(np.abs(candidate[..., column])))
if row < 2
else np.nanmax(np.abs(values))
)
image = axes[row, column].pcolormesh(
x_nm, y_nm, values.T, shading="auto", cmap="seismic", vmin=-limit, vmax=limit
)
axes[row, column].set(title=f"{row_name}: {quantity}", xlabel="x (nm)", ylabel="y (nm)", aspect="equal")
for boundary in cube.union_boundary("xy"):
axes[row, column].plot(*boundary, "k-", linewidth=1)
fig.colorbar(image, ax=axes[row, column], label=unit)
fig.suptitle(title)
return fig, axes
comparison_figure(
B_ubermag[:, :, plane_index, :],
B_microcubed[:, :, plane_index, :],
[r"$B_x$", r"$B_y$", r"$B_z$"],
f"Stray field at z={plane_z_nm:g} nm",
"T",
)
plt.show()
Gradient comparison¶
The Microcubed finite-difference result is a diagnostic, not a separate magnetic-field solver. It separates discrepancies between the underlying fields from errors introduced by approximating their derivatives on a grid.
Calculation |
Field source |
Gradient calculation |
|---|---|---|
Microcubed analytical |
Analytical cuboid field |
Analytical derivative, |
Microcubed finite difference |
Analytical cuboid field |
Numerical differences on the OOMMF grid |
OOMMF finite difference |
OOMMF demagnetizing field |
The same numerical differences |
For example, a central difference approximates
Here \(h=5\,\mathrm{nm}\). Both numerical gradients use second-order central
stencils inside the grid and second-order one-sided stencils at its outer edges.
Tensor indices are [derivative axis, field component], so [0, 2] is
\(\partial_x B_z\). All gradients below are in T/nm.
Why apply finite differences to an analytical field?¶
Comparing an analytical Microcubed gradient directly with an OOMMF numerical gradient mixes two effects: differences between the fields themselves and finite-difference approximation error. Microcubed evaluates the continuous cuboid field at observation points, whereas OOMMF returns discretized, cell-averaged fields. Neither calculation relaxes the prescribed magnetization.
Let \(D_h\) denote the common numerical differentiation operator, \(B_M\) the Microcubed field, and \(B_O\) the OOMMF field. Their gradient discrepancy obeys
The first term measures the effect of replacing the analytical derivative with a grid derivative. The second compares both fields using identical stencils and spacing. This is an exact decomposition of the differences; the individual RMS errors do not generally add, because contributions can reinforce or cancel each other. Numerical differentiation can also amplify small field errors.
If the finite-difference results agree much more closely than the analytical Microcubed and numerical OOMMF gradients, that indicates that derivative approximation at this spacing accounts for most of the larger discrepancy. It does not establish equivalence for other geometries, grids, or surface distances.
cell_nm = CELL * 1e9
dB_ubermag = np.stack(np.gradient(B_ubermag, cell_nm, axis=(0, 1, 2), edge_order=2), axis=-2)
dB_microcubed_fd = np.stack(np.gradient(B_microcubed, cell_nm, axis=(0, 1, 2), edge_order=2), axis=-2)
plane_points_nm = (
np.stack([xx[:, :, plane_index].ravel(), yy[:, :, plane_index].ravel(), zz[:, :, plane_index].ravel()]) * 1e9
)
dB_microcubed_exact = cube.dBfield(plane_points_nm).transpose(2, 0, 1).reshape(mesh.n[0], mesh.n[1], 3, 3)
dB_ubermag_plane = dB_ubermag[:, :, plane_index]
dB_microcubed_fd_plane = dB_microcubed_fd[:, :, plane_index]
assert np.isfinite(dB_microcubed_exact).all()
gradient_rows = []
for derivative, derivative_name in enumerate(AXIS_NAMES):
for component, component_name in enumerate(AXIS_NAMES):
label = f"d{derivative_name}B{component_name}"
for method, candidate in (
("microcubed analytical", dB_microcubed_exact),
("microcubed finite difference", dB_microcubed_fd_plane),
):
gradient_rows.append(
{
"component": label,
"method": method,
**error_metrics(
dB_ubermag_plane[..., derivative, component],
candidate[..., derivative, component],
),
}
)
gradient_metrics = pd.DataFrame(gradient_rows).set_index(["component", "method"])
gradient_metrics
for method in ("microcubed analytical", "microcubed finite difference"):
errors = gradient_metrics.xs(method, level="method")["NRMSE"] * 100
print(f"{method} vs. OOMMF: {errors.min():.4f}%–{errors.max():.4f}% NRMSE across nine components")
All nine analytical gradient components¶
One figure is produced for each derivative direction. As above, rows contain microcubed, Ubermag/OOMMF, and their difference.
for derivative, derivative_name in enumerate(AXIS_NAMES):
comparison_figure(
dB_ubermag_plane[..., derivative, :],
dB_microcubed_exact[..., derivative, :],
[rf"$\partial_{derivative_name} B_{component_name}$" for component_name in AXIS_NAMES],
rf"Gradient components $\partial_{derivative_name} B$ at z={plane_z_nm:g} nm",
"T/nm",
)
plt.show()
Line profiles¶
A central line on the exterior plane makes small systematic differences easier to see.
y_index = int(np.argmin(np.abs(y)))
fig, axes = plt.subplots(2, 1, figsize=(9, 7), sharex=True, constrained_layout=True)
for component, component_name in enumerate(AXIS_NAMES):
axes[0].plot(x_nm, B_microcubed[:, y_index, plane_index, component], label=f"microcubed B{component_name}")
axes[0].plot(x_nm, B_ubermag[:, y_index, plane_index, component], "--", label=f"OOMMF B{component_name}")
axes[1].plot(
x_nm, dB_microcubed_exact[:, y_index, 0, component], label=rf"microcubed $\partial_x B_{component_name}$"
)
axes[1].plot(
x_nm, dB_ubermag_plane[:, y_index, 0, component], "--", label=rf"OOMMF $\partial_x B_{component_name}$"
)
axes[0].set(ylabel="B (T)", title=f"Line at y={y[y_index] * 1e9:g} nm, z={plane_z_nm:g} nm")
axes[1].set(xlabel="x (nm)", ylabel="gradient (T/nm)")
for axis in axes:
axis.grid(alpha=0.25)
axis.legend(ncol=2, fontsize=8)
plt.show()
Interpretation¶
The tables and figures above are computed by this run, not stored benchmark claims. OOMMF returns cell-averaged demagnetizing fields; Microcubed evaluates the continuous cuboid solution at cell centers. Both use prescribed uniform magnetization, with no relaxation. Exterior fields are compared after converting H (A/m) to B (T). Global metrics exclude the magnet and a two-cell boundary layer.
The OOMMF gradient is a finite difference, whereas Microcubed also provides an analytical derivative. The additional Microcubed finite-difference result helps separate field disagreement from derivative truncation error. Gradient metrics use an exterior plane; repeat at smaller cell sizes to study convergence.
Which gradient should applications use?¶
Use dBfield() for Microcubed calculations that need gradients. It avoids
choosing a finite-difference step and subtracting nearby field values. The
finite-difference calculation here is useful for validating the analytical
derivative, comparing against grid-based solvers, and studying convergence.
In a smooth exterior region, second-order finite-difference error scales as \(O(h^2)\): halving the step should reduce that contribution by approximately a factor of four, until round-off or other errors become important. Compare at the same physical observation points and keep every stencil outside the magnet; moving closer to a surface while refining changes the problem. Near surfaces, rapid spatial variation and discontinuities can invalidate that simple trend.
A useful convergence study therefore distinguishes (1) refining the numerical derivative of the Microcubed field alone from (2) refining the OOMMF mesh, which also changes cell averaging and the underlying field calculation. The single 5 nm comparison above diagnoses agreement at one resolution; it is not itself a convergence study.