Fmm harmonic compression

benchmarks/README.md

@@ -0,0 +1,91 @@
+# Benchmarks
+
+Speed comparisons for fidimag's interactions, as opposed to the correctness
+tests under `tests/`.
+
+## `demag_fmm_vs_fft.py`
+
+Times a single `compute_field()` call for `DemagFMM` against the FFT-based
+`Demag`, on cubic `LxLxL` atomistic cuboid meshes, sweeping `L` and `theta`.
+Adapted from exploratory scripts originally written against an older
+fmmgen/fidimag integration, rerun here against the harmonic-compressed FMM
+operators (see `doc/physics_num_methods/demag_fmm.rst`).
+
+```
+.venv/bin/python3 benchmarks/demag_fmm_vs_fft.py \
+    --l-min 5 --l-max 60 --thetas 0.3 0.5 0.7 0.9 \
+    --order 8 --ncrit 128 \
+    --out benchmarks/results/demag_fmm_vs_fft.csv
+```
+
+Results are written incrementally as CSV (`L, N, order, theta, ncrit,
+t_fmm_s, t_fft_s, field_rel_err`), the last a free byproduct of timing both
+methods on the same spin configuration. Plot with:
+
+```
+.venv/bin/python3 benchmarks/plot_demag_fmm_vs_fft.py \
+    benchmarks/results/demag_fmm_vs_fft.csv \
+    benchmarks/results/demag_fmm_vs_fft.png
+```
+
+## `demag_fmm_vs_fft_sparse.py`
+
+Same comparison, but on a fixed array of nanodisks with growing spacing
+between them, rather than a growing solid cuboid. The active (`mu_s != 0`)
+site count stays fixed while the FFT method's bounding-box grid grows, which
+is the geometry the `mu_s == 0` tree exclusion in `DemagFMM` is for. See
+`doc/physics_num_methods/demag_fmm.rst`, "Where FMM/Barnes-Hut is worth it".
+
+```
+.venv/bin/python3 benchmarks/demag_fmm_vs_fft_sparse.py \
+    --n-side 4 --radius 2.5 --nz 3 \
+    --periods 6 8 10 14 18 24 32 42 55 70 90 \
+    --order 8 --thetas 0.3 0.5 0.7 0.9 --ncrit 128 \
+    --out benchmarks/results/demag_fmm_vs_fft_sparse.csv
+```
+
+Results are written incrementally as CSV (`period_nm, N_grid, N_active,
+sparsity, order, theta, ncrit, t_fmm_s, t_fft_s, field_rel_err`), the last
+restricted to the active (`mu_s != 0`) sites (see `time_demag`'s docstring
+in the script). Plot with:
+
+```
+.venv/bin/python3 benchmarks/plot_demag_fmm_vs_fft_sparse.py \
+    benchmarks/results/demag_fmm_vs_fft_sparse.csv \
+    benchmarks/results/demag_fmm_vs_fft_sparse.png
+```
+
+`plot_nanodisk_array.py` renders the `mu_s` pattern itself at a single
+period, so the geometry behind the numbers is visible:
+
+```
+.venv/bin/python3 benchmarks/plot_nanodisk_array.py \
+    --n-side 4 --radius 2.5 --nz 3 --period 18 \
+    --out benchmarks/results/nanodisk_array_geometry.png
+```
+
+## `demag_fmm_2d_vs_general.py`
+
+Times `fidimag.extensions.fmm.FMM2D` (fmmgen's planar operator variant,
+which `DemagFMM` selects automatically for a 2D, `nz=1` mesh) against the
+general `FMM` called directly on the same 2D mesh's coordinates, forcing
+the path a 3D mesh would otherwise take. See
+`doc/physics_num_methods/demag_fmm.rst`, "2D systems".
+
+```
+.venv/bin/python3 benchmarks/demag_fmm_2d_vs_general.py \
+    --l-min 10 --l-max 150 --l-step 10 \
+    --order 8 --thetas 0.3 0.5 0.7 0.9 --ncrit 128 \
+    --out benchmarks/results/demag_fmm_2d_vs_general.csv
+```
+
+Results are written incrementally as CSV (`L, N, order, theta, ncrit,
+t_planar_s, t_general_s, field_rel_err`), the last the relative L2
+difference between the two operator sets on identical input (both correct,
+so this should stay small rather than trend to zero). Plot with:
+
+```
+.venv/bin/python3 benchmarks/plot_demag_fmm_2d_vs_general.py \
+    benchmarks/results/demag_fmm_2d_vs_general.csv \
+    benchmarks/results/demag_fmm_2d_vs_general.png
+```

benchmarks/demag_fmm_2d_vs_general.py

@@ -0,0 +1,106 @@
+"""
+Benchmark: DemagFMM's automatic 2D (planar) routing vs the general 3D
+operators, on the same 2D atomistic mesh.
+
+A 2D mesh (nz=1) has every site at the same z, so fmmgen's planar operator
+variant (generate_code(..., planar=True)) can be used exactly, not as an
+approximation - it simply omits the always-zero terms the general 3D
+operators would still carry (see doc/physics_num_methods/demag_fmm.rst,
+"2D systems"). This times DemagFMM (which auto-selects the planar variant
+for a 2D mesh) against fidimag.extensions.fmm.FMM called directly on the
+same 2D mesh's coordinates, forcing the general 3D path fmmgen.extensions.
+fmm.FMM2D would otherwise replace.
+
+Usage:
+    .venv/bin/python3 benchmarks/demag_fmm_2d_vs_general.py
+"""
+import argparse
+import csv
+
+import numpy as np
+
+import fidimag
+import fidimag.extensions.fmm as fmm_ext
+import fidimag.common.constant as const
+
+A_LATTICE = 0.2715  # nm, matches demag_fmm_vs_fft.py
+MU_S = 3 * const.mu_B
+
+
+def time_one(fn, repeats):
+    import time
+    t = 0.0
+    for _ in range(repeats):
+        start = time.time()
+        fn()
+        t += (time.time() - start) / repeats
+    return t
+
+
+def time_demag(L, order, theta, ncrit, repeats=3):
+    """Build an LxL 2D atomistic mesh and time one compute_field() call for
+    the planar (FMM2D) and general 3D (FMM) operators on the same
+    coordinates. Returns (n_spins, t_planar, t_general, field_rel_err)."""
+    mesh = fidimag.common.CuboidMesh(
+        nx=L, ny=L, nz=1, dx=A_LATTICE, dy=A_LATTICE, dz=A_LATTICE,
+        periodicity=(False, False, False), unit_length=1e-9,
+    )
+    n = mesh.n
+    coords3 = (mesh.coordinates * mesh.unit_length).astype(np.float64)
+    coords2 = np.ascontiguousarray(coords3[:, :2])
+    m = np.zeros(3 * n)
+    m[2::3] = 1.0
+    mu_s = np.full(n, MU_S)
+
+    planar = fmm_ext.FMM2D(n, ncrit, theta, order, coords2, m, mu_s, 0)
+    field_planar = np.zeros(3 * n)
+    t_planar = time_one(lambda: planar.compute_field(field_planar), repeats)
+    field_planar = field_planar.copy() * -1e-7
+
+    general = fmm_ext.FMM(n, ncrit, theta, order, coords3, m, mu_s, 0, True)
+    field_general = np.zeros(3 * n)
+    t_general = time_one(lambda: general.compute_field(field_general), repeats)
+    field_general = field_general.copy() * -1e-7
+
+    field_rel_err = (np.linalg.norm(field_planar - field_general)
+                      / np.linalg.norm(field_general))
+
+    return n, t_planar, t_general, field_rel_err
+
+
+def main():
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--l-min", type=int, default=10)
+    parser.add_argument("--l-max", type=int, default=150)
+    parser.add_argument("--l-step", type=int, default=10)
+    parser.add_argument("--order", type=int, default=8)
+    parser.add_argument("--ncrit", type=int, default=128)
+    parser.add_argument("--thetas", type=float, nargs="+",
+                         default=[0.3, 0.5, 0.7, 0.9])
+    parser.add_argument("--repeats", type=int, default=3)
+    parser.add_argument(
+        "--out", type=str,
+        default="benchmarks/results/demag_fmm_2d_vs_general.csv")
+    args = parser.parse_args()
+
+    with open(args.out, "w", newline="") as f:
+        writer = csv.writer(f)
+        writer.writerow(["L", "N", "order", "theta", "ncrit",
+                          "t_planar_s", "t_general_s", "field_rel_err"])
+
+        for L in range(args.l_min, args.l_max + 1, args.l_step):
+            for theta in args.thetas:
+                n, t_planar, t_general, err = time_demag(
+                    L, args.order, theta, args.ncrit, args.repeats)
+                writer.writerow([L, n, args.order, theta, args.ncrit,
+                                  f"{t_planar:.6f}", f"{t_general:.6f}",
+                                  f"{err:.6e}"])
+                f.flush()
+                print(f"L={L:3d} N={n:6d} theta={theta}: "
+                      f"planar={t_planar:.4f}s general={t_general:.4f}s "
+                      f"(speedup={t_general / t_planar:.2f}x) "
+                      f"field_rel_err={err:.2e}")
+
+
+if __name__ == "__main__":
+    main()

benchmarks/demag_fmm_vs_fft.py

@@ -0,0 +1,107 @@
+"""
+Benchmark: DemagFMM vs the FFT-based Demag, on cuboid atomistic systems.
+
+This times a single compute_field() call for both interactions, on cubic
+LxLxL atomistic cuboid meshes of increasing size. It is a speed comparison
+only - see tests/test_demag_fmm.py for correctness checks.
+
+Adapted from exploratory scripts originally written against an older
+fmmgen/fidimag integration (test_fmm.py / test_fidimag.sh), rerun here
+against the harmonic-compressed FMM operators.
+
+Usage:
+    .venv/bin/python3 benchmarks/demag_fmm_vs_fft.py
+    .venv/bin/python3 benchmarks/demag_fmm_vs_fft.py --l-min 5 --l-max 60 \
+        --thetas 0.3 0.5 0.7 0.9 --order 5 --ncrit 128 \
+        --out benchmarks/results/demag_fmm_vs_fft.csv
+"""
+import argparse
+import csv
+import time
+
+import numpy as np
+
+import fidimag
+import fidimag.common.constant as const
+
+A_LATTICE = 0.2715  # nm, bcc Fe-like lattice constant, matches the source script
+J_EXCHANGE = 5.88 * const.meV
+MU_S = 3 * const.mu_B
+
+
+def time_demag(L, order, theta, ncrit, repeats=3):
+    """Build an LxLxL cuboid atomistic Sim and time one compute_field() call
+    for DemagFMM and for the FFT-based Demag, each averaged over `repeats`
+    calls. Returns (n_spins, t_fmm, t_fft, field_rel_err), where
+    field_rel_err is the relative L2 norm between the two methods' field --
+    a free byproduct of timing both, since both are already computed on the
+    same spin configuration."""
+    mesh = fidimag.common.CuboidMesh(
+        nx=L, ny=L, nz=L, dx=A_LATTICE, dy=A_LATTICE, dz=A_LATTICE,
+        periodicity=(False, False, False), unit_length=1e-9,
+    )
+
+    sim = fidimag.atomistic.Sim(mesh, name="demag_bench", driver="llg")
+    sim.set_mu_s(MU_S)
+    sim.add(fidimag.atomistic.Exchange(J_EXCHANGE))
+    demag_fmm = fidimag.atomistic.DemagFMM(order, ncrit, theta)
+    demag_fft = fidimag.atomistic.Demag()
+    sim.add(demag_fmm)
+    sim.add(demag_fft)
+    sim.set_m((0, 0, 1))
+
+    t_fmm = 0.0
+    for _ in range(repeats):
+        start = time.time()
+        fmm_field = demag_fmm.compute_field()
+        t_fmm += (time.time() - start) / repeats
+    fmm_field = fmm_field.copy()
+
+    t_fft = 0.0
+    for _ in range(repeats):
+        start = time.time()
+        fft_field = demag_fft.compute_field()
+        t_fft += (time.time() - start) / repeats
+    fft_field = fft_field.copy()
+
+    field_rel_err = (np.linalg.norm(fmm_field - fft_field)
+                      / np.linalg.norm(fft_field))
+
+    return mesh.n, t_fmm, t_fft, field_rel_err
+
+
+def main():
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--l-min", type=int, default=5)
+    parser.add_argument("--l-max", type=int, default=60)
+    parser.add_argument("--l-step", type=int, default=1)
+    parser.add_argument("--order", type=int, default=8)
+    parser.add_argument("--ncrit", type=int, default=128)
+    parser.add_argument("--thetas", type=float, nargs="+",
+                         default=[0.3, 0.5, 0.7, 0.9])
+    parser.add_argument("--repeats", type=int, default=3)
+    parser.add_argument("--out", type=str,
+                         default="benchmarks/results/demag_fmm_vs_fft.csv")
+    args = parser.parse_args()
+
+    with open(args.out, "w", newline="") as f:
+        writer = csv.writer(f)
+        writer.writerow(["L", "N", "order", "theta", "ncrit",
+                          "t_fmm_s", "t_fft_s", "field_rel_err"])
+
+        for L in range(args.l_min, args.l_max + 1, args.l_step):
+            for theta in args.thetas:
+                n, t_fmm, t_fft, field_rel_err = time_demag(
+                    L, args.order, theta, args.ncrit, args.repeats)
+                writer.writerow([L, n, args.order, theta, args.ncrit,
+                                  f"{t_fmm:.6f}", f"{t_fft:.6f}",
+                                  f"{field_rel_err:.6e}"])
+                f.flush()
+                print(f"L={L:3d} N={n:7d} theta={theta}: "
+                      f"fmm={t_fmm:.4f}s fft={t_fft:.4f}s "
+                      f"(fmm/fft={t_fmm / t_fft:.1f}x) "
+                      f"field_rel_err={field_rel_err:.2e}")
+
+
+if __name__ == "__main__":
+    main()