Skip to content

Fmm harmonic compression - #182

Open
rpep wants to merge 8 commits into
masterfrom
fmm-harmonic-compression
Open

Fmm harmonic compression#182
rpep wants to merge 8 commits into
masterfrom
fmm-harmonic-compression

Conversation

@rpep

@rpep rpep commented Aug 30, 2026

Copy link
Copy Markdown
Member

Pulls in new FMM perf. improvements, should be substantially faster. Also had fixed a sign bug so it'll be more accurate too.

rpep added 2 commits August 30, 2026 11:40
Ports the current fmmgen driver (Morton-ordered SoA arrays, runtime
variant selection between uncompressed and harmonic-compressed
operators) into fidimag's fmmlib, replacing a pre-S2M snapshot. The
compressed multipole/local arrays are (p+1)^2 coefficients instead of
C(p+3,3), which measures ~1.7x faster at the M2L kernel level for this
source_order=1 (point dipole) configuration.

- example.py: regenerate operators.{cpp,h} with compress=True alongside
  the existing config (order=8, source_order=1, harmonic_derivs=True).
- fmm.pyx / demag.py: DemagFMM gains a `compressed` flag (default True)
  threaded through to fmm_select(), which must run before build_tree
  since it determines the coefficient array sizes.

Two bugs in fmmgen's driver were found and fixed upstream while wiring
this up (both are pre-existing, not introduced by compression):
- Tree::compute_field_* read source values from a body_S array
  snapshotted once in build_tree, so a Tree reused across repeated
  solves with new spin/moment data each time (this integrator's usage
  pattern) kept using the values from construction. Fixed by
  re-reading through the live Particle pointers every solve.
- The parallel dual-tree traversal double-counted every P2P/M2L
  interaction on trees small/shallow enough to fully classify during
  the breadth-first expansion phase (a few hundred particles) -- the
  loop broke without swapping frontier to the now-empty next set, so
  the depth-first phase re-classified and re-appended the same pairs.
  fmmgen's own benchmarks only use N>=100k and never hit this.

Verified: DemagFMM (fmm and bh, compressed and uncompressed) agree with
brute-force DemagFull to ~1e-15 relative L2 error; full test suite
(187 passed) shows no regressions.
Existing tests only checked DemagFMM/DemagFull and Demag/DemagFull
separately, never DemagFMM against Demag directly, and only ever at
theta=0.0 (exact pairwise sum, no multipole approximation exercised).

_run_comparison now takes an optional reference interaction (defaults
to DemagFull, unchanged for existing tests). Three new tests compare
DemagFMM against Demag: 2D and 3D at theta=0.0 (near machine precision,
same as the DemagFull comparisons), plus one at theta=0.2 to actually
exercise the multipole acceptance criterion, at a tolerance matched to
the measured accuracy for that theta/order/ncrit combination.
@davidcortesortuno

davidcortesortuno commented Aug 30, 2026

Copy link
Copy Markdown
Collaborator

Looks great to me. I cannot provide precise feedback as there are lots of technical details for this massive source code. Tests look fine, which are an advantage when adding/improving new features with AI. Can you provide:

  • Benchmark testing speeds against normal FFT demag?
  • Documentation in docs (AI generated is fine, quite clever, you can follow the NEBM docs which were updated recently; I'd even pass your paper or so, so the machine summarizes what the code is doing)

I can try to add some of them if you need to

@davidcortesortuno

Copy link
Copy Markdown
Collaborator

CI failure: two bugs in fmmgen, the first hiding the second

The build fails at

fatal error: /tmp/fidimag_regen/operators.h: No such file or directory

Both causes are in fmmgen/writer.py, not in
anything hand-written here, so the fix belongs upstream and then the operators
get regenerated.

1. The include names the path the header was written to

fidimag/atomistic/fmmlib/operators.cpp:1 on this branch:

#include "/tmp/fidimag_regen/operators.h"

operators.h is committed right beside it, and master has the correct
#include "operators.h". This is the only absolute path anywhere in the
branch's build inputs.

writer.py emits #include "{name}.{hext}", and name also decides where
the files are written:

f = open(f"{src_dir.rstrip('/')}/{name}.{fext}", "w")
...
f.write(f'#include "{name}.{hext}"\n')

So generating with name='/tmp/fidimag_regen/operators' — the natural way to
put the output somewhere — writes the files correctly and poisons the include.

2. What that was hiding: __restrict at call sites

Fix the include and it gets as far as the compiler, then fails with 396
errors across 198 call sites
:

void S2M(double x, double y, double z, double * __restrict S, double * __restrict M, int order) {
switch (order) {
  case 2:
    S2M_2(x, y, z, __restrict S, __restrict M);   // not valid C++

__restrict is a declaration qualifier and cannot appear at a call site. The
dispatch is built by deleting type keywords from the declaration to turn it
into a call:

replaced_code = (
    func.replace(f"_{start}", f"_{i}")
    .replace("* ", "")        # "double * __restrict S" -> "double __restrict S"
    .replace("double ", "")   #                        -> "__restrict S"
    .replace("float ", "")
    .replace("void ", "")
)

__restrict is not in that list, so it survives. master's generated file has
zero occurrences of the pattern, because its signatures do not carry
__restrict and the string-munging happened to work.

Suggested fix

Take the last identifier of each parameter rather than deleting keywords, and
emit the basename in the include:

def _call_arguments(declaration):
    """The argument names of a C declaration, ready to be called with.

    >>> _call_arguments('void f(double x, double * __restrict M)')
    'x, M'
    """
    inside = declaration[declaration.index("(") + 1:declaration.rindex(")")]
    names = [re.split(r"[\s*]+", p.strip())[-1]
             for p in inside.split(",") if p.strip()]
    return ", ".join(names)
        called = func.split("(")[0].split()[-1]
        arguments = _call_arguments(func)
        for i in range(start, order):
            code += "  case {}:\n".format(i)
            code += "    {}({});\n    break;\n".format(
                called.replace(f"_{start}", f"_{i}"), arguments)
-    f.write(f'#include "{name}.{hext}"\n')
+    f.write(f'#include "{os.path.basename(name)}.{hext}"\n')

Verified

Regenerated with this branch's own options as recorded in operators.h
(compress=True, source_order=1, atomic=True, language='c++'), at order 4 so
it is quick, deliberately passing a directory as name to reproduce the
original bug:

line 1 of operators.cpp:   #include "operators.h"
__restrict at call sites:  0
dispatch:                  S2M_2(x, y, z, S, M);
g++ -fsyntax-only -fopenmp: compiles cleanly

I have not pushed anything to fmmgen or to this branch. Once the generator is
fixed, regenerating at order 13 should be all this PR needs — worth running in
the background, the cost climbs steeply with order.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants