Physical test scenario 6 - demagnetization field regression test to analytical cuboid solution

[AdamsMP92]

Motivation

The demagnetising field of a uniformly magnetised rectangular cuboid has a closed-form analytical solution.

This provides a quantitative reference for testing Demag.effective_field in finite-difference micromagnetic calculators.

The test would validate:

• the sign convention of the demagnetising field;
• the component ordering of the computed vector field;
• field units in (\mathrm{A,m^{-1}});
• the behaviour of Demag.effective_field for a spatially resolved sample;
• and convergence of the finite-difference demag field under mesh refinement.

Physical model

Consider an axis-aligned rectangular cuboid centred at the origin,

$$C = \left[-\frac{L_x}{2}, \frac{L_x}{2}\right] \times \left[-\frac{L_y}{2}, \frac{L_y}{2}\right] \times \left[-\frac{L_z}{2}, \frac{L_z}{2}\right],$$

with uniform magnetisation

$$\mathbf{M}_0 = \begin{bmatrix} M_x\\\ M_y\\\ M_z \end{bmatrix}.$$

The demagnetising field is

$$\mathbf{H}_\mathrm{d}(\mathbf{x}) = \frac{1}{4\pi} \Psi(\mathbf{x})\mathbf{M}_0,$$

where the geometry tensor is obtained from an alternating sum over the eight cuboid corners,

$$\Psi_{\alpha\beta}(\mathbf{x}) = \sum_{i,j,k\in\{0,1\}} (-1)^{i+j+k} \psi_{\alpha\beta} \left( \mathbf{x}; \mathbf{y}_{ijk} \right).$$

The corner positions are

$$\mathbf{y}_{ijk} = \begin{bmatrix} (-1)^i L_x/2\\\ (-1)^j L_y/2\\\ (-1)^k L_z/2 \end{bmatrix}.$$

For

$$\mathbf{r} = \mathbf{x} - \mathbf{y}_{ijk}, \qquad R = \lVert \mathbf{r} \rVert,$$

the kernel functions are

$$\psi_{\alpha\alpha} = \mathrm{atan2} \left( r_x r_y r_z, r_\alpha^2 R \right),$$

and, for (\alpha\ne\beta),

$$\psi_{\alpha\beta} = \log\left|R-r_\gamma\right|,$$

where (\gamma) is the remaining coordinate direction.

For numerical robustness, the logarithmic term can be evaluated using the algebraically equivalent expression

$$R-r_\gamma = \frac{r_\alpha^2+r_\beta^2}{R+r_\gamma}$$

in the cancellation-prone case (R\approx r_\gamma).

Proposed benchmark geometry

A small cuboid can be used, for example

lengths = (20e-9, 10e-9, 6e-9)
Ms = 1e6
M0 = (0, Ms, 0)

The corresponding Ubermag setup is

import numpy as np
import discretisedfield as df
import micromagneticmodel as mm

lengths = np.array((20e-9, 10e-9, 6e-9))
Ms = 1e6
M0 = (0, Ms, 0)

p1 = tuple(-0.5 * lengths)
p2 = tuple(0.5 * lengths)

mesh = df.Mesh(
    region=df.Region(p1=p1, p2=p2),
    n=(20, 10, 6),
)

system = mm.System(name="cuboid_demag_field")
system.energy = mm.Demag()
system.m = df.Field(
    mesh,
    nvdim=3,
    value=M0,
    norm=Ms,
)

The demag effective field can then be computed with the calculator under test,

H_demag = calculator.compute(
    system.energy.demag.effective_field,
    system,
)

or, for a system containing only the demag term,

H_demag = calculator.compute(
    system.energy.effective_field,
    system,
)

Analytical reference implementation

A compact NumPy reference implementation can evaluate the analytical cuboid field at the mesh-cell centres.

import numpy as np


def psi(alpha, beta, X, Y):
    X = np.asarray(X, dtype=float)
    Y = np.asarray(Y, dtype=float)

    r = X - Y
    R = np.linalg.norm(r, axis=-1)

    if alpha == beta:
        numerator = np.prod(r, axis=-1)
        denominator = r[..., alpha] ** 2 * R
        return np.arctan2(numerator, denominator)

    gamma = 3 - alpha - beta

    direct_argument = np.abs(R - r[..., gamma])
    transverse_square = r[..., alpha] ** 2 + r[..., beta] ** 2
    conjugate_denominator = np.abs(R + r[..., gamma])

    with np.errstate(divide="ignore", invalid="ignore"):
        stable_argument = transverse_square / conjugate_denominator

    use_stable_argument = (
        (r[..., gamma] >= 0)
        & (conjugate_denominator > 0)
    )

    argument = np.where(
        use_stable_argument,
        stable_argument,
        direct_argument,
    )

    with np.errstate(divide="ignore", invalid="ignore"):
        return np.log(argument)


def cuboid_geometry_tensor(X, lengths):
    X = np.asarray(X, dtype=float)
    lengths = np.asarray(lengths, dtype=float)

    tensor = np.zeros(X.shape[:-1] + (3, 3), dtype=float)

    for i in range(2):
        for j in range(2):
            for k in range(2):
                sign = (-1) ** (i + j + k)
                corner = np.array(
                    [
                        (-1) ** i * lengths[0] / 2,
                        (-1) ** j * lengths[1] / 2,
                        (-1) ** k * lengths[2] / 2,
                    ]
                )

                for alpha in range(3):
                    for beta in range(3):
                        tensor[..., alpha, beta] += (
                            sign * psi(alpha, beta, X, corner)
                        )

    return tensor


def cuboid_demag_field(X, lengths, M0):
    Psi = cuboid_geometry_tensor(X, lengths)
    return np.einsum("...ij,j->...i", Psi, M0) / (4 * np.pi)

The mesh-cell centres can be obtained from the mesh iterator,

points = np.asarray(list(mesh)).reshape(*mesh.n, 3)

H_reference = cuboid_demag_field(
    X=points,
    lengths=lengths,
    M0=np.asarray(M0),
)

Important interpretation

This should be a convergence or regression benchmark, not a strict pointwise equality test at a single mesh resolution.

The analytical expression gives the field of a continuously magnetised cuboid evaluated at points. A finite-difference calculator returns the demag field of a discretised, piecewise-constant magnetisation on cells. Therefore, exact pointwise agreement is not expected at finite cell size.

The expected behaviour is convergence under mesh refinement, especially away from the sample boundary.

Proposed convergence test

Run the same physical cuboid at several resolutions, for example

resolutions = [
    (10, 5, 3),
    (20, 10, 6),
    (40, 20, 12),
]

For every resolution:

1. compute Demag.effective_field;
2. evaluate the analytical cuboid field at the same cell centres;
3. compute the vector error

error = np.linalg.norm(
    H_demag.array - H_reference,
    axis=-1,
)

and the RMS error

rms_error = np.sqrt(np.mean(error**2))
relative_rms_error = rms_error / Ms

The main assertion could check that the error decreases under mesh refinement,

assert rms_error_fine < rms_error_coarse

or, preferably, that the RMS error excluding boundary cells decreases.

Boundary-cell handling

Boundary cells should either be excluded from the primary metric or reported separately.

The analytical expression is a point-value field, whereas the finite-difference demag field is associated with cell values. Close to surfaces, edges, and corners, the difference between point values, cell averages, and boundary-limit conventions is most pronounced.

A simple interior mask could exclude one cell shell from all sample boundaries,

mask = np.ones(mesh.n, dtype=bool)

mask[0, :, :] = False
mask[-1, :, :] = False
mask[:, 0, :] = False
mask[:, -1, :] = False
mask[:, :, 0] = False
mask[:, :, -1] = False

interior_rms_error = np.sqrt(
    np.mean(error[mask] ** 2)
)

The primary convergence assertion could then use interior_rms_error.

Suggested assertions

The test can verify:

1. the computed demag field has the expected shape and units;
2. all field values are finite;
3. the interior RMS error decreases under mesh refinement;
4. the field respects the expected symmetries of the cuboid;
5. selected components have the expected sign for a given magnetisation direction;
6. boundary-cell errors are reported separately from interior errors.

For magnetisation along (y),

$$\mathbf{M}_0 = M_s \mathbf{e}_y,$$

the field should obey the symmetries of the cuboid geometry. For example, in the centre plane and along symmetry lines, selected transverse components should vanish or change sign according to reflection symmetry.

Expected value

This benchmark would provide a quantitative reference for Demag.effective_field.

It would complement existing demag tests that check whether the demag term runs successfully, and would make it easier to catch regressions involving sign conventions, component ordering, field-unit handling, calculator-specific demag-field computation, and mesh-refinement behaviour.

Because the reference is analytical and independent of a second micromagnetic calculator, it could be useful both as a regression test and as a diagnostic benchmark for calculator backends.