Method of Moments for Demagnetization#
This page describes the physics and numerical method behind apply_demag,
starting from the material model and ending with all the algorithmic
optimizations that make the iterative solver fast.
1. Physical model#
A soft magnetic material responds to an applied magnetic field by acquiring a magnetization proportional to the local total field:
where \(\chi\) is the (isotropic scalar) magnetic susceptibility. In magpylib’s convention, the primary field quantity is the intrinsic polarization \(\mathbf{J} = \mu_0 \mathbf{M}\) (units of Tesla), so the constitutive relation becomes
The total field seen by a cell has three contributions:
Contribution |
Symbol |
Origin |
|---|---|---|
Remanent polarization |
\(\mathbf{J}^{(0)}\) |
Hard-magnet-like “frozen” moment |
External applied flux |
\(\mathbf{B}_\text{ext}\) |
User-specified |
Demagnetizing flux |
\(\mathbf{B}_\text{demag}\) |
The cell’s own magnetization (and neighbors’) acts back |
Self-consistency requires that the polarization satisfies this equation simultaneously for every cell:
2. Method of Moments — general formulation#
Discretization#
The magnetic body is divided into \(N\) cuboid cells, each assumed to carry a uniform polarization \(\mathbf{J}_i\) (Method of Moments / volume-element approach). Under this assumption the demagnetizing flux at cell \(i\) due to cell \(j\) is linear in \(\mathbf{J}_j\):
where \(\mathsf{T}_{ij}\) is a 3×3 demagnetization tensor block (dimensionless) that encodes the volume-averaged field produced at cell \(i\)’s centre by a unit polarization in cell \(j\).
Linear system#
Substituting into the self-consistency equation and collecting terms:
In matrix–vector form (stacking all 3-component cell vectors into a single \(3N\)-vector):
where:
\(\mathbf{J} \in \mathbb{R}^{3N}\) — unknown polarization vector (Fortran-flat: \([J_{x,1},\ldots,J_{x,N},\, J_{y,1},\ldots,J_{z,N}]\))
\(\mathbf{S} = \operatorname{diag}(\chi_1,\ldots,\chi_N,\chi_1,\ldots,\chi_N,\chi_1,\ldots,\chi_N)\) — susceptibility matrix (3N×3N diagonal)
\(\mathbf{T} \in \mathbb{R}^{3N \times 3N}\) — full demag tensor matrix (\(\mathsf{T}_{ij}\) blocks on rows \(i\), columns \(j\))
Solving this linear system gives the self-consistent polarization of every cell, accounting for all mutual demagnetization interactions.
3. The demagnetization tensor#
Newell 1993 analytical formula#
For identical, axis-aligned rectangular prisms of dimensions \((a,b,c)\), the volume-averaged interaction tensor \(\mathsf{N}_{ij}\) (a.k.a. the Newell kernel) has a closed analytical form given by
A. J. Newell, W. Williams, D. J. Dunlop (1993) A generalization of the demagnetizing tensor for nonuniform magnetization. J. Geophys. Res. 98(B6), 9551–9555.
The six independent components (\(\mathsf{N}_{xx}\), \(\mathsf{N}_{yy}\), \(\mathsf{N}_{zz}\), \(\mathsf{N}_{xy}\), \(\mathsf{N}_{xz}\), \(\mathsf{N}_{yz}\)) are computed from two auxiliary functions \(f(x,y,z)\) and \(g(x,y,z)\) via a 27-point second-difference operator:
Implemented in newell.py
(newell_f, newell_g, demag_block).
Sign convention in the code#
T_code[k, i, j, m] = -(1/mu_0) * N_mk( pos[j] - pos[i] )
After the T *= mu_0 step in apply_demag (which promotes from H-field units
to B-field units), T_code becomes dimensionless: it maps a polarization
\(\mathbf{J}\) (Tesla) to the demagnetizing flux \(\mathbf{B}_\text{demag}\)
(Tesla).
Self-demagnetization#
For the diagonal block \(i = j\) (a cell acting on itself), the Newell kernel reduces to the classical shape demagnetization factors \((N_{xx}^{(\text{self})},\, N_{yy}^{(\text{self})},\, N_{zz}^{(\text{self})})\) satisfying Brown’s identity:
For a cube: \(N_{xx} = N_{yy} = N_{zz} = \tfrac{1}{3}\).
These self-factors are used to build the Jacobi preconditioner for GMRES (see below).
4. Solvers#
4.1 Direct solver (solver="direct")#
Assemble the full \(3N \times 3N\) matrix \(\mathbf{Q} = \mathbf{I} - \mathbf{S}\,\mathbf{T}\) explicitly.
Call
numpy.linalg.solve(Q, rhs).
Cost: \(O(N^2)\) memory, \(O(N^3)\) solve. Exact to floating-point precision. Practical for \(N \lesssim 3000\).
The tensor \(\mathbf{T}\) is built via magpy.getH with three orthogonal unit
polarizations applied to all cells simultaneously (the “3 pol unit” loop in
demag_tensor).
4.2 Iterative solver (solver="iterative") — GMRES#
Instead of assembling \(\mathbf{Q}\) explicitly, GMRES only requires a matvec \(\mathbf{v} \mapsto \mathbf{Q}\,\mathbf{v}\) at each iteration. This opens the door to specialized fast matvec implementations.
The preconditioner is a diagonal (Jacobi) operator using the analytical self-demag factors:
which approximates the diagonal of \(\mathbf{Q}\) and dramatically reduces the number of GMRES iterations needed for high-susceptibility materials.
5. Optimizations#
5.1 FFT convolution for a uniform single-group grid#
When all \(N\) cells share the same dimensions and orientation, and their centres lie on a regular Cartesian grid of shape \((N_x, N_y, N_z)\) with spacing \((s_x, s_y, s_z)\), the demag tensor is translation-invariant:
The matvec \(\mathbf{H} = \mathbf{T}\,\mathbf{J}\) is then a 3-D discrete convolution, computable in \(O(N \log N)\) with an FFT on a doubled (zero-padded) grid for open boundary conditions.
Physical grid (Nx × Ny × Nz) → Doubled grid (2Nx × 2Ny × 2Nz)
_________ _________________
| | | |0 0 0 0|
| cells | | cells |0 0 0 0|
|_________| |_________|0 0 0 0|
|0 0 0 0 0|0 0 0 0|
|_________________|
The kernel (6 independent Newell components pre-FFTed) is built once by
build_fft_kernel and reused for every GMRES iteration.
Flow for a single-group FFT matvec:
J (3N flat, original-cell order)
│
▼ reorder by grid permutation
J_grid (Nx, Ny, Nz, 3)
│
▼ zero-pad to doubled grid
J_pad (2Nx, 2Ny, 2Nz, 3)
│
▼ rfftn
J_fft (2Nx, 2Ny, Nz+1, 3) [half-spectrum]
│
▼ multiply by kernel FFT (Nxx,Nyy,Nzz,Nxy,Nxz,Nyz)
H_fft (2Nx, 2Ny, Nz+1, 3)
│
▼ irfftn → crop to (Nx, Ny, Nz, 3)
│
▼ reorder back to original-cell order
H (3N flat)
│
▼ return J - S * H [the (I - ST) matvec result]
Implemented in demag_fft_matvec and _build_fft_matvec.
5.2 Block-FFT for multi-group assemblies#
When the input is a collection of several meshed cuboids (each with a
different orientation or size), no single uniform grid spans all cells. Instead,
detect_grid_groups partitions the cells into groups:
Group: a maximal subset of cells sharing the same cell dimensions and the same orientation (quaternion), whose centres form a valid regular Cartesian grid.
In practice, each meshed cuboid produces exactly one FFT group.
The full \(3N \times 3N\) demag matrix is then naturally split into a block structure:
Groups: G₁ G₂ G₃
┌─────────┬─────────┬─────────┐
G₁ │ T₁₁ FFT │ T₁₂ CSR │ T₁₃ CSR │
├─────────┼─────────┼─────────┤
G₂ │ T₂₁ CSR │ T₂₂ FFT │ T₂₃ CSR │
├─────────┼─────────┼─────────┤
G₃ │ T₃₁ CSR │ T₃₂ CSR │ T₃₃ FFT │
└─────────┴─────────┴─────────┘
Diagonal blocks T_GG → handled by group-local FFT kernel
Off-diagonal blocks T_AB → handled by sparse CSR matrices (see §5.3)
The block matvec for one GMRES step is:
For each group G:
H_G += FFT_self_block(J_G) ← group-local FFT kernel
for each other group A ≠ G:
H_G += T_AG (CSR) @ J_A ← sparse cross-block matvec
Rotation handling: each group may have an arbitrary orientation \(R_G\). The FFT operates in the group’s local frame, so the polarization is rotated into that frame before the FFT and the result is rotated back:
J_global → R_G⁻¹ → J_local → FFT → H_local → R_G → H_global
Implemented in _fft_block_h (per-group FFT self-block) and
_build_block_fft_matvec (full block matvec).
5.3 Sparse triage for cross-block interactions#
For two groups \(A\) and \(B\) separated by a large distance, most cell pairs \((i \in A,\, j \in B)\) have negligible mutual influence. The \(\chi \cdot V / r^3\) triage criterion exploits this:
where
\(\varepsilon = \) solver_tol * 0.1 and \(\mathbf{r}_{ij} = \mathbf{p}_i - \mathbf{p}_j\)
is the cell-centre displacement.
Physical interpretation: \(V_j / r^3\) is proportional to the solid angle subtended by cell \(j\) as seen from cell \(i\). The \(\chi\) weighting accounts for the fact that a low-susceptibility cell scarcely responds even to a strong incident field.
The critical distance below which no pairs are dropped:
For 128-element 1 mm³ cubes with \(\chi = 1\) and solver_tol = 1e-6:
\(r_\text{triage} \approx 200\,\text{mm}\).
Algorithm in _compute_cross_block:
Input: n_a cells in group A, n_b cells in group B
1. Compute influence(i,j) = max(χᵢ,χⱼ) * Vⱼ / |rᵢⱼ|³ for all (i,j) [O(n_a·n_b)]
2. Build boolean mask = influence > ε
3. Prune active sources:
active_b = {j : mask[:,j].any()} ← columns with ≥1 significant pair
active_a = {i : mask[i,:].any()} ← rows with ≥1 significant pair
4. Call magpy.getH only for (active_a × active_b) sub-problem [reduced cost]
5. Zero out remaining below-threshold entries within the sub-block
6. Scatter into full CSR matrix of shape (3·n_a, 3·n_b)
The resulting T_AB is stored as a scipy.sparse.csr_matrix. Each GMRES matvec
then costs \(O(\text{nnz})\) instead of \(O(n_a \cdot n_b)\).
6. Complete algorithm flow#
apply_demag(collection, solver="iterative")
│
├─ 1. Collect cells: positions, dimensions, orientations, χ
│
├─ 2. Detect grid structure
│ ├─ detect_uniform_grid → single FFT group?
│ │ YES → fft_info path
│ ├─ detect_grid_groups → multiple FFT groups?
│ │ YES → block-FFT + sparse cross-blocks path
│ └─ otherwise → dense T matrix (magpy.getH)
│
├─ 3. Build demag operator
│ ├─ [single group] build_fft_kernel → kernel_fft stored in fft_info
│ ├─ [multi-group] build_fft_kernel per group
│ │ _compute_cross_block for every (A,B) pair → CSR T_AB
│ └─ [dense] demag_tensor via magpy.getH → T (3N×3N dense)
│
├─ 4. Solve linear system (I - S T) J = J⁽⁰⁾ + S B_ext
│ ├─ [direct] numpy.linalg.solve(Q, rhs) O(N³)
│ └─ [iterative] GMRES with Jacobi preconditioner
│ matvec options:
│ · _build_fft_matvec (single group) O(N log N)
│ · _build_block_fft_matvec (multi-group) O(N log N + nnz)
│ · _build_dense_matvec (fallback) O(N²)
│
└─ 5. Assign self-consistent polarizations back to source objects
7. Complexity summary#
Configuration |
Tensor build |
Matvec per GMRES iteration |
|---|---|---|
Single uniform grid |
\(O(N \log N)\) |
\(O(N \log N)\) |
\(K\) groups, well-separated |
\(O(N \log N) + O(\text{nnz})\) |
\(O(N \log N) + O(\text{nnz})\) |
\(K\) groups, touching |
\(O(K^2 N^2 / K^2) = O(N^2)\) |
\(O(N^2)\) |
Dense fallback |
\(O(N^2)\) |
\(O(N^2)\) |
where nnz \(\ll K^2 (N/K)^2\) when groups are separated beyond \(r_\text{triage}\).
8. References#
Newell, A. J., Williams, W., & Dunlop, D. J. (1993). A generalization of the demagnetizing tensor for nonuniform magnetization. Journal of Geophysical Research: Solid Earth, 98(B6), 9551–9555.
Chadbec, O. et al. (2006). Micromagnetic simulation using the demagnetization tensor approach. (Method of Moments references in
demag_tensordocstring.)