We need answer user: "Design a 2d Nbody gravity simulation to run on an nvidia gpu. Use the fast multiple method. Return an explanation of your approach and python pseudo code." Need produce final with explanation and Python pseudo code. Need likely correct method is Fast Multipole Method (FMM) maybe user typo "fast multiple method". Need design 2D N-body gravity on NVIDIA GPU, using FMM. Provide approach: spatial decomposition tree/box, multipole expansions for far fields, local expansions / particle-particle near fields, GPU parallelism: one thread per particle, shared memory for small cell, CUDA graphs? Python pseudo code using PyTorch CUDA or CuPy maybe pseudocode. Need maybe include equations, expansions: Newtonian potential Phi = G sum m_j/|r_i - r_j|. In 2D? Usually N-body in 2D physical plane with 3D Newtonian gravity? Or 2D gravity? "2d Nbody gravity simulation" likely x,y positions but inverse distance in 2D? Gravity in 2D could be ln? In Newtonian gravity restricted to 2D plane with masses in plane uses 3D r = sqrt(dx^2+dy^2), potential -G/r. If truly 2D gravitational potential, force ~ 1/r, potential ~ ln r. Need mention choose kernel; FMM standard for 1/r. Could specify 2D coordinates, 3D kernel; if you want 2D log kernel, replace expansions. Need likely answer robust. Need design: - Inputs N particles positions x,y, masses, dt, theta accuracy, max degree L. - Build quadtree each step? FMM needs tree. For adaptive GPU: build quadtree on CPU/GPU? For pseudo code, use uniform grid or static quadtree? To keep GPU efficient, use uniform grid of cells / tree? FMM can be implemented with uniform grid and multi-level (M3?) But easiest: regular octree/quadtree with max depth, one thread per leaf? However dynamic tree per timestep costly. Could use fixed quadtree domain; if particles move, assign to leaves. Need near interactions: within same or neighbor cells (27? for 2D quadtree 9 neighboring cells). For far fields, compute cell multipole moments and local expansions. Need pseudo code in Python: maybe using cupy or torch. Provide pseudocode with annotations. Need include M2L conversion for 2D. Maybe simpler: compute far force directly using multipole sum: For target particle i in cell B, for well-separated source cell A, contribution approximated by expanding kernel in Taylor about centers: 1/|r_i - r_j| = sum_{l=0}^L (1/l!) (D_i^l)(1/R) * ( (r_j - c_A)^l )? More precisely multipole moment M_A,l = sum_j (r_j - c_A)^l m_j; far force on target at r: F = -G sum_A sum_l M_A,l ∇_r (1/|r - c_A|) derivatives evaluated at? Need sign. Potential Phi_i ≈ G sum_A sum_l M_A,l ∇_r^l (1/R)?? Let's derive. Kernel K(r_i, r_j) = 1/|r_i - r_j|. Let R = r_i - c_A, δ = r_j - c_A. 1/|R - δ| = sum_{|l|<=L} (δ^l / l!) (∂^l / ∂R^l) 1/|R|? Taylor f(R-δ)=sum (-δ·∇)^n f(R)? Let's test 1D: 1/(x-δ)=1/x * 1/(1-δ/x)=1/x sum δ^n/x^n = sum δ^n/x^{n+1}. Derivatives ∂^n (1/x) = (-1)^n n! / x^{n+1}. So ∑ δ^n/n! *? If use ∂^n: ∂^n f = (-1)^n n!/x^{n+1}; to get positive need (-1)^n δ^n /n! ∂^n f = δ^n/x^{n+1}. So expansion: 1/|R-δ| = sum_l (-1)^l/l! (δ^l) (∂^l_R 1/R). Or equivalently sum (δ^l/l!) (∂^l_{-R} 1/R)? In FMM standard multipole: Φ(x) ≈ Σ M_l ∇^l 1/|x - c|. Maybe signs not essential in pseudo code if implemented with standard formulas. For 2D, can use complex coordinate z=x+iy. Expansion: 1/|z - ζ| = Re(1/(z-ζ))? Actually 1/sqrt((x-ξ)^2+(y-η)^2). Complex expansion for 1/(z-ζ) not same; kernel is real part of 1/(z-ζ)? 1/(z-ζ) = (x-ξ + i(y-η))/r^2, real part = (x-ξ)/r^2, not 1/r. So not simple. Could implement Cartesian tensor derivatives. Pseudocode can abstract `M2L` or `M2P` kernels. But user asks design, explanation, python pseudo code. Need not exact formula? Could provide formulas. Need GPU design details: - Use one thread per particle to compute forces; memory coalesced; shared memory for cell moments; vectorized; double precision optional. - FMM tree: 1. Domain bounding box, max depth D, leaf size max particles or min cell size. Build quadtree. For dynamic, use occupancy counting and scan; or use fixed tree. 2. For each node, compute multipole moments up to order p (e.g. 4-8) from children or particles. `M[l] = Σ m (r - c)^l`. 3. Traversal: For each particle, traverse tree; if node well-separated (not overlapping target cell and distance > θ * sizes), use M2P evaluate multipole; else recurse; if leaf near, do direct P2P. Alternative bottom-up FMM: compute local expansions via M2L for far cells, then L2P for particles; near interactions in leaves. More efficient O(N log N) or O(N). For GPU, bottom-up uniform grid simpler: For each level, compute multipole moments per cell; M2L to local cells; L2P; near P2P. But quadtree traversal per particle may be okay. Need maybe include Barnes-Hut? "fast multiple method" maybe Fast Multipole Method. We can say assume FMM (Fast Multipole Method); if you meant Fast Multipole/Barnes-Hut. Use multipole expansions. Need Python pseudo code. Could present two pseudocode styles: - High-level `class GPUNBodyFMM` with `build_tree`, `compute_multipole`, `m2p`, `step`. - Use `cupy` or `torch` pseudo code. Need maybe "python pseudo code" not fully working. Could use `@cuda.jit`? Simpler: use PyTorch CUDA with vectorized operations for fixed grid. But dynamic tree in torch hard. Pseudocode can be conceptual. Maybe include a compact pseudocode with a GPU kernel `force_kernel` per particle traverses a static quadtree stored in arrays; plus a setup kernel to assign particles to leaves. Need ensure it's plausible. Let's design: Assumptions: - 2D positions (x,y), masses, Newtonian 1/r kernel (in the plane). If true 2D logarithmic, swap kernel derivatives. - Use quadtree over square domain [0,L]^2 with max depth `D`. Each node has up to 4 children. Node id integer; children 4*id+1..+4? Need tree arrays. - Use fixed max depth; if a leaf overflows, increase depth? For pseudo code can use `max_depth` enough; or rebuild. - Multipole degree `p` e.g. 4. Moments: M[node, l1, l2] for |l1|+|l2|<=p? Number K = (p+1)(p+2)/2. For p=4 K=15. For p=6 K=28. Use double for accuracy. - Direct near interactions for particles in same or neighboring leaves (8-neighborhood) or nodes whose boxes are not well-separated. Algorithm options: A. Particle-based traversal (top-down) per particle: - For particle i, starting root, if node contains i (target cell) direct interactions with all particles in node? Need avoid double counting? If each target particle traverses tree, for source nodes not overlapping target cell and well separated use M2P; if near, descend to leaves and P2P with particles in leaves. But if we do for every target particle and every source leaf, we double count? Force on i due to j computed; j due to i computed separately; that's okay (Newton's third not exploited). Need skip self. But if near interactions for leaf with target i, we sum over all particles j in leaf except i. If we traverse for every i, each i-j near pair computed twice (i->j and j->i) but that's just separate target forces, not duplicate in same force? Actually force on i due to j computed once in i's traversal; force on j due to i computed in j's traversal. That's fine. No double count for same target. However if node near target cell, we descend and sum over all particles in near source leaves. For each target i, near pairs computed O(n_near). Good. But if a target particle in leaf, source leaf same, we sum all in leaf. If target i, source leaf near, we sum all particles in leaf. For every target, near cost maybe O(N * avg near particles). If leaf size bounded, okay. - Need well-separated criterion: node box does not intersect target cell expanded? Use standard: if node box is outside an annulus around target cell: distance from target center to source center > θ * (size_source + size_target?) More precise for quadtree: source node well-separated from target cell if source box does not touch target cell's adjacent cells? For Barnes-Hut: well separated if source node is not one of 3x3 neighbor cells of target cell. In bottom-up, use 9 neighbor cells per level. In top-down, if target leaf known, source node is well-separated if its box is not overlapping target leaf or adjacent leaves? Could use: distance between centers > θ * (size_source + size_target) or `not adjacent`. Simpler: Use θ criterion: `dist(c_i, c_node) > θ * (width_i + width_node)`. But if target in node? skip. This is approximate. - M2P evaluation: For source node A with center c_A, radius w_A, multipole moments M_A up to p. For target position r, R = r - c_A. Compute potential or force via derivatives. Pseudo code can call `evaluate_m2p(M, R, w_A, p)`. Need formulas for force. We can provide generic tensor derivative: Potential approximation: Φ_i ≈ G Σ_A Σ_{|l|≤p} M_{A,l} (-1)^{|l|} / l! (∂^l_R (1/|R|))_{R=r_i-c_A} Force f_i = -∇_r Φ_i = -G Σ_A ...? Let's derive sign. Actual potential Φ(r)=G Σ m/|r - r_j| (positive? gravitational potential often -Gm/r, force = -m ∇Φ? If Φ=-Gm/r, force = -m ∇Φ = -Gm r/r^3? Let's choose force formula direct: Force on i due to j: F_ij = G m_i m_j (r_j - r_i) / |r_j - r_i|^3 (attractive). In code can omit m_i and store acceleration a_i = G Σ m_j (r_j - r_i)/r^3. For multipole, approximate acceleration: a_i ≈ G Σ_A Σ_l M_{A,l} * ? derivative of kernel w.r.t r_i? Actual a_i = G Σ_j m_j ∇_{r_j} (1/|r_i-r_j|)? Because ∇_{r_i} (1/r) = -(r_j-r_i)/r^3? Let's check: K=1/|r_i-r_j|. ∇_{r_i} K = -(r_i - r_j)/r^3 = (r_j-r_i)/r^3. Yes. So a_i = G Σ m_j ∇_{r_i} K. If expansion K ≈ Σ M_l C_l(R), with R=r_i-c_A, M_l = Σ m_j (-1)^l? Then a_i = G Σ M_l ∇_r C_l(r_i-c_A). We can define `M2P_accel`. Could provide formula for derivative of 1/R: ∇ (1/R) = -R/R^3? Wait R=r_i-c_A. If source at c_A, K ≈ 1/|R|. ∇_r 1/|R| = -R/|R|^3. But actual force is G m (r_j-r_i)/r^3 = G m (-(R - δ))/|...|^3. If approximate source at center: a ≈ -G M0 R/R^3, but should point to source center: R = r_i - c, -R/R^3 points to c. Good. Taylor for acceleration: K(r_i,r_j)=1/|R-δ| with δ = r_j-c. M_l maybe Σ (-δ)^l? Let's set: 1/|R-δ| = Σ_{l} (1/l!) (δ·∇_δ?)? Let's choose convention: M_A,l = Σ_j (-1)^{|l|} m_j δ^l / l! ? Then K ≈ Σ_l M_A,l (∂^l_R 1/R)? Let's test 1D. Want K=1/(R-δ)=Σ δ^n/R^{n+1}. Derivatives ∂^n_R (1/R)=(-1)^n n!/R^{n+1}. If define M_n = (-1)^n δ^n/n! * m? Then M_n ∂^n = (-1)^n δ^n/n! * (-1)^n n!/R^{n+1}=δ^n/R^{n+1}. So M_n = Σ m (-1)^n δ^n/n!. Then K=Σ M_n ∂^n (1/R). Then acceleration a=G Σ ∇ K ≈ G Σ M_l ∇ ∂^l_R (1/R). In 1D: ∇ ∂^n 1/R = derivative w.r.t R: ∂^{n+1} 1/R. For n=0: ∂ 1/R = -1/R^2. a = G M0 ∂ 1/R = -G m / R^2. If R positive target right of source, acceleration negative (left), good. For n=1: M1 = -m δ; ∂^2 1/R=2/R^3; product -2mδ/R^3, expansion of -1/R^2? derivative actual -1/(R-δ)^2 ≈ -1/R^2 - 2δ/R^3. yes. So acceleration = G Σ_A Σ_l M_{A,l} ∂^{l+1?} Actually vector derivative of ∂^l scalar: ∇ (∂^l 1/R). Can compute. But in pseudocode can abstract. Maybe easier: M2P computes acceleration using recursive formulas for derivatives of inverse distance. Provide function: ``` def inv_r_derivatives(R, p): # returns D[l] where D[l] = ∂^l (1/|R|) ``` Then `a += G * M[node] * grad_D[l]`. But need not implement all. Could also use local expansion to reduce repeated M2P for many particles in same cell: L2P. But top-down M2P per particle per well-separated node may be O(N log N). Bottom-up local expands can be O(N). The user specifically asks FMM; local expansions are core. We should include M2L and L2P. Maybe design with uniform grid for simplicity and GPU friendly. Option: Multi-level uniform grid FMM (also called FMM). Steps: - Create grid at levels l=0..L, cell size h_l = 2^{L-l} maybe. - For each level: - Assign particles to finest cells; propagate to coarser cells? Or compute multipole moments for each cell directly by summing particles in cell (O(N log L) with multi-level). - Compute M2L: for each cell B, for each well-separated cell A at same level (not in 3x3 neighborhood), transfer A's multipole to B's local expansion. Use translation T_{l}. - Compute L2P: for each particle in cell B, evaluate local expansion to get acceleration. - Near field: P2P for particles within neighboring cells (3x3) at same level; possibly also across levels? In standard FMM, near interactions within leaf cells and their neighbors at same finest level? For multi-level uniform grid, if near cells are at different levels? There is a standard tree-based FMM. Uniform grid with one level can be O(N^2) if many particles in near cells. Multi-level uniform grid uses cells at multiple levels; near interactions between particles in neighboring cells at all levels? Might be complex. Simpler: static quadtree leaf level; near only same/neighbor leaves; far via M2L among nodes at same depth? If all leaves same depth (fixed tree), near 9 cells. Far cells well separated. This is essentially uniform grid if all cells at max depth. For N large, cell size small; if particles distributed, near okay, but if clustered, leaf overflow. Could allow variable depth but M2L across variable depths complex. For pseudo code, a static quadtree with variable depth and top-down M2P may be acceptable and easier. But to truly use FMM, need local expansions? Top-down M2P still multipole but not local expansion; sometimes called fast multipole? It avoids P2M per target? It's multipole expansion, but FMM includes M2L. The phrase "fast multiple method" maybe could mean "fast multipole method". Need include M2L to be safe. Let's design a GPU-friendly FMM using fixed quadtree at depth D (uniform leaf grid) with local expansions: - Domain decomposed into C x C square grid of cells at leaf level (D log2 C). Each cell has fixed capacity `Pmax`. If too many particles, increase C or use subcells; pseudo code can assume capacity. - For each cell, compute multipole moments from particles in cell (P2M). Number cells C^2. GPU: one thread per cell or block per cell; loop particles in cell (using CSR lists). Or one thread per particle atomicAdd to cell moments? Atomics can be heavy but okay; better sort particles by cell then block reduce. - M2L: For each target cell b, loop over all source cells a that are well-separated (not in 3x3 neighborhood). Naively C^4 per level; for C large, too much. Standard FMM M2L uses tree/neighbor lists and translates only well-separated cells; number O(C^2 log C) or O(C^2). For uniform grid, each cell has O(C^2) well-separated cells, total O(C^4)? Wait number cell pairs well-separated is O(C^4)? Actually number cells O(C^2), each has O(C^2) far cells, so O(C^4). But standard FMM uses multi-level to reduce. At a single level, M2P for each target cell and far source cell is O(C^4) which for C~sqrt(N) is O(N^2). Need multi-level or tree. In standard FMM, multi-level tree reduces to O(N). For pseudo code, can still mention use multi-level quadtree and M2L. But code simpler if top-down per particle O(N log N) with M2P, no local expansion. But maybe FMM expected O(N log N) acceptable? "Fast multiple method" maybe not strict O(N). Need present robust approach. Could design adaptive quadtree with local expansions using standard bottom-up traversal. But implementing in pseudo code is harder. We can outline high-level and provide simplified pseudo code for tree-based FMM. Let's recall standard quadtree FMM: - Each node has center c, size h. - Multipole moments M, local expansion L, particle list if leaf. - Well-separated: node A and B are well-separated if A not in 3x3 neighborhood of B (including diagonal) at same or different depth? Standard: if neither is ancestor of other and their boxes do not touch? In tree FMM, for each cell, far nodes are those well-separated; M2L can be done via traversal. - Algorithm: 1. P2M: bottom-up compute M for all nodes from children/particles. 2. M2L: top-down? For each node B, for each node A well-separated from B, translate M_A to L_B. Efficient traversal: for each B, traverse tree; if node A is well-separated from B, add M2L(A,B); else if A is ancestor? skip; else recurse. This is O(N log N) maybe. 3. L2P: top-down: For each particle in leaf B, add L_B evaluation (M2P of local expansion) and near P2P with particles in neighboring leaves. But to be truly fast, M2L via list of well-separated nodes can be generated with tree traversal; for pseudo code okay. GPU parallelism: - Build tree: maybe on CPU initially, update every k steps; or GPU using occupancy sort. For pseudo code, assume tree arrays known and updated. - P2M: parallel over nodes: each block computes moments for node by reducing particles in subtree? Could use one thread per leaf particle atomicAdd to ancestors? For depth D small, each particle contributes to D ancestors; atomicAdd to moments. Or bottom-up: leaf moments from particles, then internal sum children. GPU: parallel over internal nodes, 4 children vectorized. - M2L: parallel over target nodes; each block loops over well-separated source nodes from precomputed CSR `far_cells[target]`. Translation kernels: for each source-target pair, compute translation matrices T_{l}^{m} and L_B += T M_A. Number of far pairs O(N log N). Can use shared memory for source moments, vectorized. If too many, use multi-level local expansions. - L2P: parallel over particles; each particle loops over its cell's local expansion (size K) to add acceleration; plus direct near interactions: each block handles a leaf and neighbor leaves using shared memory; thread per particle in block; loop over neighbor particles. Avoid self/double count by computing all pairs or using Newton third? Simpler compute all target particles against near source leaves, skip self. If block per leaf, for target particle i in leaf B, source leaves = B and 8 neighbors. Sum. To avoid double counting within same leaf? If computing acceleration for all particles in B due to all in B, each i loops over all j in B; that's fine (force on each i). Complexity O(P_leaf^2). If using symmetric pair within block, can compute half and atomic, but not necessary for pseudo. - Integrator: Velocity Verlet: compute a(t), update x += v dt + 0.5 a dt^2, compute new a, update v += 0.5(a_old+a_new)dt. Or leapfrog. Need mention energy drift, softening epsilon to avoid singularities: K = 1/sqrt(r^2 + eps^2), derivative modified. Accuracy parameters: - multipole order p (e.g. 4 or 6) - opening angle θ (e.g. 0.5 or 1.0) controls well-separated and error. For 1/r kernel, error O((δ/R)^{p+1}). Choose p and θ such that (θ/(1+θ))^{p+1} small. Maybe mention. - softening ε. Data structures: - Particles: x, y, vx, vy, ax, ay, m; float64 or float32. On GPU arrays. - Tree node arrays: center_x, center_y, half_size, left/top child ids, first_particle, last_particle, parent, node_level, is_leaf. - Moments: M[node, K], L[node, K]; K=(p+1)(p+2)/2. - Near neighbor list: for each leaf, 9 neighboring leaf ids. - Far list: CSR of well-separated source nodes for each target node for M2L, or use traversal kernel. Need Python pseudo code. Could use PyTorch or CuPy. Maybe pseudo code with `torch.cuda` and `torch.autograd`? Pseudocode not actual. Use `@cuda.jit` style? Let's design a readable pseudo code: ``` import torch from torch import tensor # assume tensors on cuda class NBodyFMM: def __init__(self, N, domain, max_depth, p=4, theta=0.7, eps=1e-3): ... def build_quadtree(self): ... def assign_particles_to_leaf(self): ... def p2m(self): ... def m2l(self): ... def l2p(self): ... def near_p2p(self): ... def compute_accelerations(self): ... def step(self, dt): ... ``` But need more detailed pseudocode kernels. Maybe present "Python pseudo-code (not optimized, conceptually runnable)". Use a simplified static uniform leaf grid for clarity, then note for full FMM replace with multi-level quadtree. But user asked FMM, so maybe use quadtree with local expansion. Let's combine: explanation with full standard FMM, pseudo code with simplified quadtree far list. Potential issue: "fast multiple method" maybe "fast multipole method" but maybe "Fast Multipole Method" in 2D. Need maybe say "I assume you mean Fast Multipole Method (FMM)". Need final answer: explanation of approach and Python pseudo code. Should be thorough. Let's formulate mathematical expansions. 2D multipole moments: For node A with center c_A, particle j position r_j, δ_j = r_j - c_A. For multi-index α=(α1, α2), |α|=α1+α2≤p: M_{A,α} = Σ_{j∈A} m_j (-1)^{|α|} δ_j^{α} / α! where α! = α1!α2!, δ^α = δ_x^{α1} δ_y^{α2}. Potential/acceleration far evaluation (M2P): R = r_i - c_A. Define S(R)=1/|R| (or softened). For each α, compute G_{A,α}(R)=∂^α S(R) and its gradient ∇ ∂^α S(R). Then acceleration contribution: a_i^{far} += G m_i? Actually acceleration = G Σ m_j (r_j-r_i)/r^3. If M includes m_j, then a_i += G Σ_α M_{A,α} ∇ ∂^α S(R) (with sign convention above). Need verify sign with α multi-index. In 1D, M_n = (-1)^n m δ^n/n!, S=1/R? Actually |R| positive; ∂^n |R|^{-1} for R>0 = (-1)^n n! / R^{n+1}. M_n ∂^n = δ^n/R^{n+1}. Good. For vector R, derivative includes sign based on components; standard Taylor: f(R-δ) = Σ (δ^α/α!) (∂^α_{-R} f)(R) = Σ (-1)^{|α|} δ^α/α! (∂^α_R f)(R). So M_α = Σ m (-1)^{|α|} δ^α/α!. Then f=Σ M_α ∂^α_R f. Acceleration = G Σ ∇_r f = G Σ M_α ∇ ∂^α_R f. Good. For 0 order: M0=m, ∇ S = -R/R^3; acceleration = -G m R/R^3 = G m (c - r)/R^3, points to source center. Good. Local expansion (L): For target cell B with center c_B, we want approximate forces on particles i∈B from far source A. Let η = r_i - c_B. Source positions relative to A center? M2L converts M_A to L_B such that: 1/|r_i - r_j| ≈ Σ_β L_{B,β} ∂^β_{?} K? Standard local expansion: f(r_i) ≈ Σ_β L_{B,β} ∂^β_{c_B} f? Let's derive. For target near c_B, R = r_i - c_A = (c_B - c_A) + η. Let D = c_B - c_A. Expand f(R) in η around D: f(D+η) = Σ η^β/β! (∂^β_D f)(D). But f(D+η) = K(r_i, source). If M_A approximates source: K ≈ Σ_α M_{A,α} ∂^α_R S(R). Need combine: K ≈ Σ_α M_{A,α} ∂^α_R S(D+η) = Σ_α M_{A,α} Σ_β η^β/β! (∂^{α+β}_D S)(D). So L_{B,β} = Σ_{α: |α|+|β|≤p?} M_{A,α} ∂^{α+β}_D S(D). But local expansion can be evaluated as Σ_β L_{B,β} η^β? For acceleration: a_i = G Σ ∇_r K ≈ G Σ_β L_{B,β} ∇ η^β = G Σ_{β: |β|>0} L_{B,β} e^{β}? If L defined for potential. But if we want force directly, define local force coefficients. Simpler: M2L computes local potential coefficients L_B, then L2P computes potential and gradient (force). Use formulas. Local evaluation: Φ_i^{far} ≈ G Σ_B Σ_β L_{B,β} (r_i - c_B)^β? If L_{B,β} = Σ_A Σ_α M_{A,α} ∂^{α+β}_{D} S(D). Then acceleration = -? If potential positive? We'll just say L2P evaluates local expansion and its gradient using same derivative routines. Maybe avoid exact sign by saying "using standard FMM translation kernels". For pseudo code, define: ``` def M2P_accel(M, R): # M: moments of source cell, R=target - source center return G * sum_{alpha} M[alpha] * grad_partial_S(R, alpha) def M2L_coeffs(M_A, D, p): # D = c_B - c_A return sum_{alpha,beta} M_A[alpha] * partial_S(D, alpha+beta) def L2P_accel(L, eta): # eta = r_i - c_B return G * sum_{beta} L[beta] * grad_eta_beta # or gradient of local potential ``` Need be careful with signs and p split: For multipole degree p, in M2L we can use degree p. If local degree p. For each target cell B, sum over source A: for α,β with |α|+|β|≤p, L_B,β += M_A,α * ∂^{α+β} S(c_B-c_A). But if M_A only up to p, local coefficients up to p. In L2P, use up to p. Error O((cell size/target distance)^{p+1}). But for force directly, if L is potential coefficient, force = -∇ potential? In gravitational potential often Φ = -GΣm/r; force = -∇Φ? Let's avoid by defining `acceleration` kernels directly. Could define source multipole coefficients for acceleration: M^F_{α} = Σ m (-1)^{|α|} δ^α/α! M2P acceleration: G Σ M^F_α ∇ ∂^α S. M2L for acceleration: For target cell B, local coefficients C_β such that a_i ≈ G Σ_β C_β ∇ η^β. Derive: a(r_i)=G Σ_α M_α ∇ ∂^α S(D+η). Expand ∂^α S(D+η)=Σ_β η^β/β! ∂^{α+β} S(D). Then a=G Σ_β [Σ_α M_α ∂^{α+β} S(D)] ∇ η^β /? Wait expansion includes η^β/β!, and ∇ η^β = β_γ η^{β-eγ}/(β-eγ)!? Better store local coefficients including 1/β!. Let C_β = Σ_α M_α (∂^{α+β} S(D))/β!. Then a=G Σ_β C_β ∇(η^β). For β=0, ∇(1)=0 so no contribution (constant potential). Good. M2L: ``` for alpha, beta with |alpha|+|beta| <= p: L_B[beta] += M_A[alpha] * partial_S(D, alpha+beta) / factorial(beta) ``` L2P: ``` a += G * sum_{beta, |beta|<=p, beta != 0} L_B[beta] * grad(monomial_eta(beta)) ``` This works. M2P can also be viewed local with eta? M2P acceleration: G Σ_α M_α ∇ ∂^α S(R). Need derivative routine for softened inverse distance in 2D. For S(R)=1/sqrt(R_x^2+R_y^2+eps^2). Need partial derivatives up to order p. Could compute via recurrence: Let r2 = R_x^2+R_y^2+eps^2, r = sqrt(r2). For first derivatives: ∂_x S = -R_x / r^3. ∂_y S = -R_y / r^3. Second: ∂_xx S = (3R_x^2 - r^2)/r^5? Check: derivative -x r^{-3}: -r^{-3} + 3x^2 r^{-5}; with r^2=x^2+y^2+eps^2 => (-r^2+3x^2)/r^5. ∂xy = 3xy/r^5. Higher via recursion: ∂_x^{a+1}∂_y^b S = -R_x ∂_x^{a}∂_y^b (1/r)? Actually derivatives of r^{-1}. Recurrence: D^{α+e_x} S = -R_x D^α S - r^2 D^{α+2e_x} S? Not simple. Use generating function or symbolic precomputed formulas. In pseudo code: `derivative_S(R, alpha, p)` can be precomputed polynomial over powers of r. For p small (4-6) hardcode or use recurrence: For f = r^{-1}, r^2 = x^2+y^2+eps^2. Identity: (x ∂_x + y ∂_y) r^{-1} = -r^{-1}. Also Euler. General recurrence: ∂_x^{a+1} ∂_y^b f = -x ∂_x^a ∂_y^b f - y^2? Hmm. Known recurrence for Cartesian derivatives of 1/r: ∂_x^{l+1}∂_y^m f = -x ∂_x^l∂_y^m f - r^2 ∂_x^{l+2}∂_y^m f? No. Let's derive: ∂_x (r^{-1}) = -x r^{-3}. More generally, for q = r^{-1}, ∂_x q = -x q / r^2. But ∂_x^{a+1} ∂_y^b q = ∂_x_y^b (q) = ∂_y^b (∂_x q) = ∂_y^b(-x q/r^2) = -x ∂_y^b(q/r^2) because x independent y. Not simple. Alternative use complex? For 2D maybe use polar/harmonic expansions. But pseudo code can say use a precomputed table or automatic differentiation. Maybe use expansion in terms of distance and angle? In 2D, kernel 1/r can be expanded in multipoles using solid harmonics. For 2D, use 2D solid harmonics (cos nθ, sin nθ). Multipole moments per order n: For source cell A: M_{A,n}^{cos} = Σ_j m_j r_j^n cos(n φ_j) M_{A,n}^{sin} = Σ_j m_j r_j^n sin(n φ_j) Far field potential in target local polar (ρ,θ): 1/R = Σ_{n=0}^∞ (r^n/R^{n+1}) cos(n(φ-θ))? Need if target cell? Actually for source cluster: if target far from cell, use source multipole: potential at target position relative to source center: 1/|R-δ| ≈ Σ_{n=0}^L (δ^n/R^{n+1}) cos(n(φ_δ - φ_R))? Let's derive using Legendre/Gegenbauer? In 2D, expansion of 1/|r-r'| in polar for r>r': 1/|r-r'| = (1/r) Σ_{n=0}∞ (r'/r)^n P_n(cos(θ-θ'))? Actually 2D uses cos nΔθ: 1/√(r^2+r'^2-2rr'cosΔ) = (1/r) Σ_{n=0}∞ (r'/r)^n cos(nΔ). Yes. This is simpler. Multipole coefficients: M_n^c = Σ m r_j^n cos(n φ_j), M_n^s = Σ m r_j^n sin(n φ_j). Then potential at target (R,Θ): Φ≈ G Σ_n (M_n^c cos nΘ + M_n^s sin nΘ)/R^{n+1}. Force = -∇Φ? If potential positive? We can compute acceleration directly: a = G ∇ (1/R)?? Let's use this for explanation. It is 2D multipole expansion. For target local expansion, expand source far field around target cell center using similar expansion: local coefficients from translation. This is more elegant for 2D. Need choose 2D kernel 1/r (in plane). The expansion is valid for 2D Euclidean distance (not 3D z=0? It is same as 3D distance in 2D). Good. Let's derive force from potential coefficients. If potential `Phi = G Σ m / R` (positive), gravitational acceleration for a test mass? Usually force = -G m_i m_j r_vec / r^3. If potential U = -G m_i m_j / r, force = -∇U = -G m_i m_j? Let's define `acceleration = G Σ m_j (r_j-r_i)/r^3`. This equals `G ∇_{r_i} (1/|r_i-r_j|)` as above. So if we compute `K = 1/R`, acceleration = G ∇_r K. If approximate K = Σ_n A_n(R,Θ), acceleration = G ∇ K. In polar, K ≈ Σ_{n=0}^L (A_n cos nΘ + B_n sin nΘ) R^{-(n+1)}. Gradient in polar: ∂_R K = -Σ (n+1)(A cos + B sin) R^{-(n+2)} (1/R)∂_Θ K = -Σ n(-A sin + B cos) R^{-(n+2)} Convert to x,y: ∇K = (cosΘ ∂_R - sinΘ/R ∂_Θ, sinΘ ∂_R + cosΘ/R ∂_Θ). Can compute. For pseudo code, can use Cartesian derivative routine. Maybe use solid harmonics for M2L translation formulas in 2D. Translation from multipole to local: Local expansion for target cell B: K ≈ Σ_{n=0}^L (C_n cos nθ + D_n sin nθ) ρ^n, where ρ,θ relative to target center. For far source cell A, given multipole moments M_n^c,s relative to A center, distance d=c_B-c_A with polar (D,Φ), translation: C_n = D^{-(n+1)} Σ_{k=0}^n M_k? Need formula. Use addition theorem. Multipole source: K ≈ Σ_k M_k^c,s (cos/sin kφ_source?) / |R|^{k+1} where R = r - c_A. For target near c_B, R = D + η (vector from A to target point? Actually R = r_i - c_A = (c_B - c_A) + η = D + η, with |η| small, |D| large). Expand 1/|D+η| in powers of η for |η|<|D|. In 2D: 1/|D+η| = (1/D) Σ_{m=0}∞ (ρ/D)^m cos(m(θ-Φ)) maybe with sign? Since |D+η|^2 = D^2+ρ^2+2Dρ cos(θ-Φ). For ρr': (1/r) Σ (r'/r)^n cos n(θ-θ'). If r=D, r'=η, angle difference θ_η - θ_D. Distance |D - η|? Our vector D+η = D - (-η). If r'= -η, angle θ' = θ+π, cos n(θ'-θ_D)=cos n(θ-Φ+π)=(-1)^n cos n(θ-Φ). So 1/|D+η| = (1/D) Σ (-ρ/D)^n cos n(θ-Φ). Yes. So local coefficients include (-1)^n. Then combine with source multipole expansion: Source term for order k: M_k^c cos kφ + M_k^s sin kφ over |D+η|^{k+1}. Need expansion of (D+η)^{-(k+1)} as local series. In 2D, (1/|D+η|^{k+1}) = D^{-(k+1)} Σ_{m=0}∞ (ρ/D)^m P? Not just cos m; for powers, use Gegenbauer polynomials. In 2D, expansion: 1/|D+η|^{k+1} = D^{-(k+1)} Σ_{m=0}∞ (ρ/D)^m C_m? Maybe `P_m^{(1/2)?}`? Actually generating function for Gegenbauer: 1/√(1 - 2t z + t^2) = Σ C_n^λ(t) z^n. For 1/(1 - 2t cosΔ + t^2)^{(k+1)/2} = Σ t^m P_m^{(k)}(cosΔ) (Gegenbauer with λ=(k-1)/2? In 2D? Let's see. For 1D? Maybe too complex. Cartesian tensors easier to describe; no need exact M2L formula in final, just standard translation kernels. But perhaps include enough. Alternative: Use "fast multipole method" in GPU with PyTorch and `scipy`? No, need pseudo. Need decide final structure: 1. Assumptions and kernel. 2. High-level design: quadtree, FMM stages, GPU mapping. 3. Mathematical model: moments, M2P/M2L/L2P, near P2P, softening, integrator. 4. GPU performance design: memory layout, parallel kernels, shared memory, atomics, double/float, tree rebuild, occupancy, profiling. 5. Python pseudo code. 6. Complexity/accuracy notes. Need ensure pseudo code is Python-like. Could use `# pseudo-code` comments. Include class with methods and kernel-like loops. Use `torch` pseudo code for GPU arrays. For FMM far list, maybe precompute `far_nodes[target]`. Let's craft pseudo code: ``` # Python pseudo-code (conceptual; use CuPy/PyTorch/CuTe for real CUDA) import torch class NBody2D_FMM: def __init__(self, pos, mass, domain, max_depth, p=4, theta=0.7, eps=1e-3, leaf_max=64): self.x = torch.tensor(pos[:,0], dtype=torch.float64, device='cuda') ... self.tree = QuadTree2D(domain, max_depth, leaf_max) self.tree.assign(self.x, self.y) self.K = (p+1)*(p+2)//2 self.M = torch.zeros(num_nodes, K, dtype=torch.float64, device='cuda') self.L = torch.zeros(num_nodes, K, dtype=torch.float64, device='cuda') self.far_lists = self.tree.build_far_lists(theta, max_degree=p) self.near_lists = self.tree.build_near_leaf_lists() def p2m(self): # parallel over leaves; then bottom-up leaf_block = ... for each leaf cell b: M[b, :] = 0 for i in leaf.particles: delta = pos[i] - center[b] for alpha <= p: M[b, idx(alpha)] += mass[i] * (-1)^|alpha| * delta^alpha / alpha! for level from max_depth-1 downto 0: for node n: M[n] = sum M[child] # child centers? Need moments about node center, not just sum child moments about child centers. ``` Important: If child moments are about child centers, parent moments about parent center require translation. In bottom-up P2M, can compute leaf moments about leaf center, then parent moments about parent center need M2C translation (multipole center shift). Simpler: compute moments directly from particles for every node? O(N log N) but okay. Or for internal nodes, use children with translation. Pseudo code can ignore translation? Need be correct. We can say use `P2M` by direct summation in each node's particle range if tree stores particles sorted by subtree; but internal node's particle range includes all descendants; one block per node loops over its particles and computes moments about node center. Complexity O(N log N) maybe but for pseudo okay. Better: use bottom-up with center translation: If child moment M_c,α about child center c_c, parent center c_p. δ_parent = (c_c - c_p) + δ_child. M_p,β = Σ_{α+γ=β} M_c,α (c_c - c_p)^γ / γ! maybe with sign? Because M_c,α = Σ m (-1)^{|α|} δ_child^α/α!. Need expand. But can just state "using center-shift translation" or compute directly. For pseudo, simplify: `M[n] = sum over particles in subtree` maybe with GPU loop. But performance? We can note real implementation uses bottom-up translation. Maybe data layout: tree is a complete quadtree; each node has particle list. For P2M, use one thread per particle atomicAdd to all ancestor moments. Since depth D small (e.g. 10), each particle does D*K operations. This is GPU-friendly and avoids child center translation? It still needs delta to each ancestor center; compute. Atomics to double can be slow; use per-block reductions. Pseudo: ``` for each particle i: node = leaf[i] while node != -1: delta = pos[i] - center[node] for alpha: M[node, alpha] atomic_add(mass[i]*...) node = parent[node] ``` But if each particle contributes to every ancestor, internal M correct. Complexity O(N D K). Good for pseudo. But FMM usually O(N log N). We can present bottom-up as preferred, but pseudo can use direct ancestor accumulation for clarity. However GPU atomic double heavy; mention in production use sorted particles and block reductions. M2L stage: ``` def m2l(self): L.zero_() # kernel: one block per target node for each target node b in parallel: for a in far_lists[b]: D = center[b] - center[a] if well_separated(a,b): # translation matrix T(D,p) for beta in 0..p: for alpha in 0..p-alpha: L[b, beta] += M[a, alpha] * partial_S(D, alpha+beta) / factorial(beta) ``` Need if local expansion for target node from source node a. But if target node has children, local expansion for parent also needed for particles in descendants? Standard FMM: local expansions at all levels; L2P passes local expansion down: L2L from parent to child. In simplified pseudo, if we compute L for every node from far lists, particles in leaves can use L of leaf and ancestors? If M2L computes L for all nodes, far fields from nodes well-separated to each target node. A particle in leaf receives contributions from far sources to leaf's L, and also from sources far from ancestor but not leaf? Actually if source is far from ancestor, it's also far from leaf? The standard hierarchical approach uses local expansions at each level to avoid missing contributions. If we only use leaf L from all far source nodes, far source list for leaf may include many nodes; still correct if all far sources not in near list are included. Could be many. For performance, use M2L at all levels and L2L. Pseudo can mention standard: compute L for every node; then propagate local expansions from parent to child (L2L) so descendants inherit far-field from ancestors. But if we already sum all far source nodes to leaf, not needed. For full FMM, include L2L. Let's outline standard bottom-up FMM stages: 1. P2M: compute M for all nodes. 2. M2L: for each node, sum M2L from well-separated nodes (not ancestors/descendants? Actually source nodes well-separated from target node). This produces L for each node from far cells at same level? In hierarchical, M2L is done node-to-node for well-separated pairs; then L2L propagates local expansion from parent to children, because a particle in child receives far fields from nodes that are well-separated from parent but maybe not child? Need standard algorithms. 3. L2P: evaluate L at particles and near P2P. We can state: For each leaf, total far field = L_leaf plus contributions from L of ancestors propagated (or use L2L). Pseudo can simplify: after M2L, do `l2l` top-down: `L[child] += translate_local_to_local(L[parent], center_child - center_parent)`. Then L2P uses leaf L. Need near P2P: ``` def near_p2p(self): ax.zero_(); ay.zero_() # one block per leaf for leaf b: neighbors = near_leaves[b] # 9 cells for each particle i in b: for each nb in neighbors: for each particle j in nb: if i==j: continue dx = x[j]-x[i]; dy=y[j]-y[i] r2=dx*dx+dy*dy+eps*eps inv_r3 = r2**(-1.5) ax[i] += G*m[j]*dx*inv_r3 ay[i] += G*m[j]*dy*inv_r3 ``` This double counts? For acceleration of i, okay. If leaf b and nb both in neighbors, when block b computes i in b due to j in nb; block nb computes j due to i. That's separate. No double count for same i because for a given i, source nb included once. If neighbors list includes b and all 8. Good. Complexity O(sum leaf neighbor pairs). If leaf max P, okay. Could optimize by block loads neighbor particles in shared memory. L2P: ``` def l2p(self): for each particle i in leaf b: eta = pos[i] - center[b] for beta in 0..p, beta != 0: ax[i] += G * L[b,beta] * d_eta_x(monoid(eta,beta)) ay[i] += G * L[b,beta] * d_eta_y(monoid(eta,beta)) ``` But local expansion coefficients might be for acceleration, as defined. Need mention if using potential, take negative gradient depending sign. Far M2P alternative: If no local expansion, for each particle traverse far nodes. But we use L. Need include M2L translation pseudo code maybe with generic `partial_S`. We can define `monomial` and `derivative_S` as helper functions. For p small, precompute translation matrices per relative offset? In uniform grid, relative center distances are discrete; precompute T for each relative offset and order to speed. For quadtree, variable distances; compute on fly or precompute per level/relative position. Could include a simplified Cartesian derivative function: ``` def partial_S(R, alpha): # returns ∂^alpha (1/sqrt(Rx^2+Ry^2+eps^2)) # use recursive formula or precomputed symbolic expressions ``` Pseudo code can leave as black box. But maybe user wants approach, not full implementation. Need mention 2D true logarithmic gravity if desired: - If "2D gravity" means infinite line masses / 2D Poisson, kernel is log r, force ∝ (r_j-r_i)/r^2. Replace S(r)=ln r and derivatives; FMM still works with log kernel (2D FMM uses logarithmic expansions). Provide note. Need maybe include "fast multiple method" could be "Fast Multipole Method" and "multiple" typo. State assumption. Let's think about GPU details deeper: - Memory coalescing: particle arrays SoA (structure-of-arrays). Tree nodes SoA. Moments 2D array (node, K) coalesced if node major. For M2L, far_lists CSR: each target node has list of source nodes; block loops over list, loads M[a] (K doubles) coalesced across threads if each thread handles one β/α? Use vectorized `float64` maybe less parallelism; use float32 for performance if tolerance. Use `torch.float64` for accuracy; can use mixed precision with error correction. - Shared memory: In near P2P, load positions/masses of leaf+neighbors into shared memory, thread per particle; if Pmax small (32-128). In M2P/L2P, one thread per particle loads L[leaf] into shared/register; K small. - Tree building on GPU: Use occupancy counting: for max depth D, compute leaf index for each particle by rounding position to leaf grid; sort by leaf using radix sort; count particles; if > leaf_max, rebuild with deeper tree or use subcells. For adaptive quadtree, GPU can use prefix sums to assign node ids. But rebuilding every step can be expensive; do every 10-100 steps or when max cell occupancy exceeds threshold. Or use fixed domain grid with periodic boundaries? Mention. - Parallel reduction: P2M: sort particles by leaf; use block per leaf to compute moments; then bottom-up internal nodes. If using ancestor atomicAdd, use `torch.add`? In CUDA use `atomicAdd` for float; double atomics less efficient; use warp shuffles and block reductions. - FMM far list: Build once per tree; store in CSR. Number far pairs can be large; if too large, use multi-level M2L and L2L to keep O(N log N). For pseudo code, far list per node. - OpenMP/CUDA graphs: capture kernels for repeated steps; use pinned memory; avoid host sync. Use `cupy` or PyTorch `torch.compile`? Not necessary. Need maybe include code with PyTorch pseudo kernels? Could write actual vectorized PyTorch for near direct? But FMM kernels in PyTorch hard. Pseudo code with loops is fine. Let's design pseudo code more concretely: ``` # ---- helpers (pseudo) ---- def monomial(delta, alpha): return delta.x**alpha[0] * delta.y**alpha[1] def factorial(alpha): return math.factorial(alpha[0]) * math.factorial(alpha[1]) def multi_index_list(p): return [(a,b) for a in range(p+1) for b in range(p+1-a)] # S(R) = 1/sqrt(Rx^2+Ry^2+eps^2) # derivative_S(R, alpha) returns partial^alpha S # grad_derivative_S(R, alpha) returns (d/dx, d/dy) partial^alpha S class QuadTree2D: def assign(self, x, y): ... def build_near_lists(self): ... def build_far_lists(self, theta): # for each node b, list nodes a such that boxes are well separated: # not ancestor/descendant and dist(center_a, center_b) > theta*(size_a+size_b) # optionally only same level? for simplified code. ``` But if far lists include all well-separated nodes for each node, M2L can be O(N^2) in worst? But tree well-separated pairs O(N log N)? For quadtree, number of well-separated node pairs is O(N log N) maybe. Let's check: In a quadtree with N points, number of well-separated pairs O(N log N) or O(N)? Standard FMM O(N). For pseudo okay. M2L code: ``` def m2l(self): self.L.zero_() # CUDA kernel pseudo: for b in range(num_nodes): # parallel for a in far_lists[b]: D = self.center[b] - self.center[a] # only if a and b are well separated (prechecked) T = translation_matrix_M2L(D, self.p, self.eps) # (K_beta x K_alpha) self.L[b, :] += T @ self.M[a, :] # L2L top-down: propagate local expansions to children for level in range(0, max_depth): for node in nodes_at_level(level): # parallel for child in children[node]: d = center[child] - center[node] self.L[child, :] += L2L(self.L[node, :], d, self.p, self.eps) ``` Need M2L translation matrix shape: L_beta = Σ_alpha T_{β,α} M_alpha. T computed from derivative_S(D, α+β)/β! maybe. L2L translation: If parent local expansion L_parent for potential/acceleration, child local L_child += translate local coefficients. For acceleration coefficients C, translation maybe C_child_γ = Σ_β L_parent_β ∂^β (monomial?)? If local expansion is polynomial in η around parent center: a(η_parent) ≈ Σ_β C_β ∇ η_parent^β. For child center offset d, η_parent = d + η_child. Then C_child_γ = Σ_β C_β ∂^β (d+η_child)^? /? This is a polynomial translation. Pseudo can call `L2L`. But for pseudo, maybe skip L2L and say in simplified version we evaluate M2P for each particle from far nodes. However user asked FMM; local expansion better. Maybe present two pseudo code variants: "Simplified GPU FMM using per-particle multipole evaluation (easier to understand)" and "Production bottom-up FMM stages". But final might be too long. Need enough. Could write pseudo code for bottom-up with black-box translation functions; that's acceptable. Let's define indices and K: ``` ALPHAS = [(a,b) for a in range(p+1) for b in range(p+1-a)] IDX = {alpha: i} K = len(ALPHAS) ``` Moments: ``` def p2m_kernel(): for i in range(N): node = leaf_of[i] while node >= 0: d = (x[i]-cx[node], y[i]-cy[node]) for a in range(K): alpha=ALPHAS[a] M[node,a] += mass[i] * ((-1)**(alpha[0]+alpha[1])) * d[0]**alpha[0]*d[1]**alpha[1] / fac[alpha] node = parent[node] ``` But doing while loop per particle on GPU in CUDA is okay if depth small. Pseudo uses atomic add; in real use block reductions. Could mark `atomic_add`. M2P evaluation (if using local? maybe not). We can include `m2p` for completeness: ``` def m2p_accel(node, r): R = r - center[node] a = 0 for alpha: dS = derivative_S(R, alpha) grad_dS = grad_derivative_S(R, alpha) a += M[node, alpha] * grad_dS return G*a ``` But with local expansion, L2P similar. Local expansion evaluation: If L is acceleration coefficients: `L[β]` such that a = G Σ L[β] ∇ η^β. Then: ``` def l2p_accel(node, r): eta = r - center[node] a=0 for beta in ALPHAS: if beta == (0,0): continue # grad of eta^beta: gx = beta[0] * eta[0]**(beta[0]-1) * eta[1]**beta[1] if beta[0]>0 else 0 gy = beta[1] * eta[0]**beta[0] * eta[1]**(beta[1]-1) if beta[1]>0 else 0 a += L[node,beta] * (gx, gy) return G*a ``` But if L coefficients include 1/β! as defined, gradient of η^β is correct. For β=(1,0), ∇x=1; L_{1,0} * 1. Does this match? Let's test 1D first-order local. Source at D>0, target near 0. K≈1/(D+η)=1/D - η/D^2 + η^2/D^3... Acceleration G ∂_η K = G(-1/D^2 + 2η/D^3 - ...). Our C_0=0? C_1 = M0 ∂^{1} S(D)/1! ? M0=m. ∂ S(D) = -1/D^2. L1 = -m/D^2. L2P: β=1: G * L1 * ∂η^1 = G * (-m/D^2). Good. β=2: M0 ∂^2 S(D)/2! = (2/D^3)/2 = 1/D^3. Gradient of η^2 = 2η. G*(1/D^3)*2η = 2Gη/D^3. Good. So works. M2L formula for acceleration coefficients: ``` L_beta += M_alpha * partial_S(D, alpha+beta) / factorial(beta) ``` where D = c_B - c_A. Check β=1, α=0: partial_S(D,1) = -1/D^2; L1 = m*(-1/D^2)/1! good. β=2: partial_S(D,2)=2/D^3; /2! =1/D^3 good. M2P formula for acceleration: `a = G Σ_alpha M_alpha * grad(partial_S(R, alpha))`. Check α=1: M1 = -mδ. partial_S(R,1) = ∂_R? In 1D ∂ S = -1/R^2. Gradient = ∂_x = 2/R^3. product -2mδ/R^3, good. In 2D, grad of partial derivative. Good. Need derivative of S with multi-index in 2D. Could implement using recursive polynomial. Let's see if we can provide a simple recurrence for Cartesian derivatives of `r^{-1}` up to p. It would make pseudo code more concrete. Let's attempt to derive a recurrence that is easy: Let f = r^{-1}, r^2 = x^2+y^2+ε^2. We can compute all derivatives D_{a,b} = ∂_x^a ∂_y^b f up to p using recurrence from differentiating `x ∂_x f + y ∂_y f = -f`? That gives relations but not full. Another recurrence: For f = (x^2+y^2+ε^2)^{-1/2}, derivatives can be expressed as polynomial in x,y times r^{-(2n+1)}. We can compute via automatic differentiation on symbolic polynomials. In pseudo code, we can say precompute `D[a,b]` using recurrence: ``` # Let r2 = x^2+y^2+eps^2, r = sqrt(r2) # For a+b > 0: # D[a,b] = -x * D[a-1,b] - y^2 * D[a-2,b]? Not correct. ``` Let's derive using `∂_x f = -x r^{-3} = -x f / r^2`. If we know D[a,b] = f * P_{a,b}(x,y) where P polynomial odd degree a+b? For f=r^{-1}. Then ∂_x (f P) = f_x P + f P_x = -x/r^2 f P + f P_x = f (P_x - x P/r^2). But r^2 = x^2+y^2+ε^2. To maintain polynomial P times f? P_x - xP/r^2 has denominator r^2; not polynomial times f? But maybe can express with r^{-3}. Better express D[a,b] = r^{-(a+b+1)} Q_{a,b}(x,y) where Q polynomial degree a+b. Recurrence for Q. Let D[a,b] = r^{-(a+b+1)} Q_{a,b}. Base Q[0,0]=1. For a>=1: ∂_x r^{-(a+b+1)} = -(a+b+1) x r^{-(a+b+3)} ∂_x Q_{a-1,b} (degree a+b-1) D[a,b] = ∂_x (r^{-(a+b-1)} Q_{a-1,b})? Wait D[a-1,b] = r^{-(a+b-1)} Q_{a-1,b}. Differentiate: D[a,b] = -(a+b-1) x r^{-(a+b+1)} Q_{a-1,b} + r^{-(a+b-1)} ∂_x Q_{a-1,b}. Factor r^{-(a+b+1)}: Q_{a,b} = ∂_x Q_{a-1,b} * r^2 - (a+b-1) x Q_{a-1,b}. Since r^2 = x^2+y^2+ε^2, Q remains polynomial degree (a+b). Good. Similarly for b>=1: Q_{a,b} = ∂_y Q_{a,b-1} * r2 - (a+b-1) y Q_{a,b-1}. This recurrence works? Test a=1,b=0: Q[1,0] = ∂_x Q[0,0]*r2 - 0*x*1 = r2 = x^2+y^2+ε^2. D[1,0]=r^{-3} r2? That would be 1/r, not -x/r^3. Missing negative sign? Because derivative of r^{-(a+b+1)}? Let's recalc: D[0,0]=r^{-1}. For a=1: D[1,0]=∂_x r^{-1} = -x r^{-3}. Formula using D[0,0] = r^{-(0+0+1)} Q[0,0]. Differentiate: derivative r^{-1} = -1 x r^{-3}. General: D[a,b] = ∂_x D[a-1,b] with D[a-1,b]=r^{-(a+b-1)} Q. Derivative: -(a+b-1) x r^{-(a+b+1)} Q + r^{-(a+b-1)} ∂_x Q. Factor r^{-(a+b+1)}: Q - (a+b-1)x Q. For a=1: Q = 1 -0 =1, D = r^{-3}, wrong sign. Because exponent for D[a-1,b] is (a-1+b+1)=a+b, not a+b-1? Let's set n=a+b. D[a,b] = ∂^n r^{-1} should be r^{-(n+1)} Q degree n. Base n=0: r^{-1}. If D[a-1,b] has n-1 derivatives: exponent = -(n). So D[a-1,b]=r^{-n} Q. Differentiate: -n x r^{-n-2} Q + r^{-n} ∂Q = r^{-(n+1)} (∂Q*r^2 - n x Q). For a=1,b=0: n=1? Wait D[0,0] n=0, derivative D[1,0] n=1. Formula using previous n=0: exponent previous = -1? Actually D[0,0]=r^{-(0+1)}=r^{-1}. Derivative: -1*x*r^{-3}. General previous n=0: D_prev = r^{-(n+1)} Q, derivative: -(n+1)x r^{-(n+3)} Q + r^{-(n+1)} ∂Q = r^{-(n+3)}(∂Q*r^2 - (n+1)x Q). New exponent should -(n+3)? New derivatives n+1, exponent -(n+2). Inconsistent because derivative of r^{-(n+1)} gives r^{-(n+3)}, but new D should r^{-(n+2)}. Wait derivative of r^{-1} should r^{-3}, but formula for n=0 says new exponent -3, but expected n+1=1 derivative exponent -3? General derivative of r^{-1} first derivative is r^{-3}, not r^{-2}. For second derivative exponent -5. So exponent = -(n+2*? Actually each derivative increases power by 2? For 1D, d^n r^{-1} = (-1)^n (2n-1)!! r^{-(2n+1)}. So each derivative increases exponent by 2. For multi-index total n=a+b, exponent = -(2n+1). My earlier assumption r^{-(n+1)} wrong. So Q polynomial degree n? First derivative D=x? D = -x r^{-3}; exponent -3, Q=-x degree1. Good. Recurrence: D[a,b] = r^{-(2n+1)} Q_{a,b}, n=a+b. Base Q[0,0]=1, D=r^{-1}. If D[a-1,b] has n-1 derivatives: exponent = -(2(n-1)+1)=-(2n-1). Derivative: ∂_x (r^{-(2n-1)} Q) = -(2n-1) x r^{-(2n+1)} Q + r^{-(2n-1)} ∂Q = r^{-(2n+1)} ( ∂Q * r^2 - (2n-1) x Q ). Thus Q[a,b] = r2*∂_x Q[a-1,b] - (2(a+b-1)+1) x Q[a-1,b] = r2 ∂_x Q[a-1,b] - (2(a+b)-1) x Q[a-1,b]. For a=1,b=0: Q = r2*0 - (1)*x*1 = -x. D = -x r^{-3}, correct. Second x: a=2,b=0, n=2, previous Q[1,0]=-x, ∂x=-1. Q = r2*(-1) - (3)x(-x)= -r2+3x^2 = 2x^2 - y^2 - ε^2. D = r^{-5}(2x^2 - y^2 - ε^2). Compare known ∂xx 1/r = (3x^2 - r^2)/r^5 = (2x^2-y^2-ε^2)/r^5, correct. Great. So we can provide recurrence for Q polynomials. Need gradient of D[a,b] for M2P/L2P? We can compute `grad_derivative_S(R, alpha)` by differentiating Q recurrence? Or simply note `grad ∂^α S = (∂_x^{α_x+1}∂_y^{α_y} S, ∂_x^{α_x}∂_y^{α_y+1} S)`. Ah! The gradient of ∂^α S is just the derivatives with α+e_x and α+e_y. So we don't need separate; compute all derivatives up to p+1 maybe. For acceleration M2P: `grad_partial_S(R, α) = (D[α+e_x], D[α+e_y])` where D[a,b]=∂_x^a∂_y^b S. So need derivatives up to p+1. For L2P local polynomial gradient is algebraic. Thus define `derivative_S(R, a,b)` using Q recurrence up to `p+1`. This is nice. Let's write helper pseudo code: ``` def make_S_derivatives(p): # returns function D(R, a, b) # Q[a,b] polynomial? To evaluate at R, we can compute recursively numerical values Qval[a,b] def D(R, a, b): n = a+b r2 = Rx*Rx + Ry*Ry + eps*eps # compute Qval[a,b] recursively memoized if a>0: Q = r2 * dQ_dx(a-1,b) - (2*n - 1)*Rx*Qval(a-1,b) elif b>0: Q = r2 * dQ_dy(a,b-1) - (2*n - 1)*Ry*Qval(a,b-1) else: Q=1 return Q / r2**(n+0.5) # because exponent 2n+1 => r2^(n+0.5) ``` But need dQ_dx, dQ_dy. We can avoid derivative of Q by using finite? For pseudo, say use symbolic polynomial table. Or use recurrence for D directly? We can derive direct recurrence for D values without Q derivatives? Q derivative needed. Could compute Q polynomials symbolically offline for p<=6. In pseudo code, `partial_S` can be a precomputed table of formulas. That's acceptable. Maybe simpler: Use automatic differentiation or precomputed formulas in real code. Pseudo code can state `# implemented with a small precomputed polynomial table or AD`. But we can provide a direct recurrence for D values using derivative of Q? Let's see if D recurrence possible: D[a,b] = ∂_x D[a-1,b]. If D[a-1,b] = r^{-(2n-1)} Q. Need Q. Could compute Q and its derivatives symbolically. For p small, fine. Alternative use polar expansion to avoid high-order Cartesian derivatives. But M2L formulas also complex. Cartesian black-box okay. Maybe provide pseudo code with `kernel_derivs` as placeholder: ``` def inv_r_derivatives(R, max_order): """ Return dict D[a,b] = ∂_x^a ∂_y^b (1/sqrt(Rx^2+Ry^2+eps^2)) for 0 <= a+b <= max_order. In production, unroll for p=4/6 or use CuPy polynomial table. """ ``` Then use it. Need M2L translation matrix with derivative up to p. For each β, α, need `D[α+β]` up to total p. So max_order p. For M2P acceleration, need gradient of D[α], i.e. D[α+e] up to p+1. So compute up to p+1. Let's define `accel_from_m2p`: ``` def m2p_contrib(M, R): D = inv_r_derivatives(R, p+1) ax=ay=0 for alpha in ALPHAS: ax += M[alpha] * D[alpha[0]+1, alpha[1]] ay += M[alpha] * D[alpha[0], alpha[1]+1] return G*(ax,ay) ``` Check sign: D[1,0] = ∂_x S = -x/r^3. M0=m => ax = G*m*(-x/r^3), where x = target - source. Points to source if x positive, negative. Good. For α=(1,0): M1 = Σ m (-δ_x). D[2,0] = ∂xx S = (3x^2-r^2)/r^5. Contribution to ax: M1*D[2,0]. In 1D, M1=-mδ, D[2]=2/R^3? With softening? If source at 0, R=r, D[2]=2/r^3, product -2mδ/r^3. Actual acceleration -Gm/(r-δ)^3? Expansion -Gm/r^3(1+3δ/r+...)? Wait earlier using K expansion: K=1/(r-δ)=1/r + δ/r^2 + δ^2/r^3; acceleration ∂K/∂r = -1/r^2 - 2δ/r^3 - 3δ^2/r^4. But multipole M1=-mδ, D[2]=∂^2(1/r)=2/r^3 => -2mδ/r^3. Good. D[2] in 1D = 2/r^3. In 2D formula D[2,0] at y=0 = (3r^2-r^2)/r^5=2/r^3. Good. M2L translation for acceleration coefficients: ``` def m2l_contrib(M_A, D): # returns L_beta S = inv_r_derivatives(D, p) L = zeros(K) for beta in ALPHAS: for alpha in ALPHAS where |alpha|+|beta| <= p: L[beta] += M_A[alpha] * S[alpha+beta] / factorial(beta) return L ``` Where `D = c_B - c_A`. Check β=0: L0 = M0*S[0,0]=m/D. But β=0 local coefficient corresponds constant acceleration? L2P ignores β=0. However local potential has constant; no force. Fine. We could store only β !=0. If using acceleration coefficients, L0 unused. But M2L to child local expansions may need β=0? L2L of acceleration coefficients maybe β=0 doesn't contribute to force but may contribute after translation to child? Constant acceleration field has zero? Wait constant acceleration (zero order local) could exist? In our local acceleration expansion, β=0 term would be `L0 * ∇1 = 0`, so it never contributes. But when translating local expansion from parent to child, a constant potential term (β=0) in parent local expansion could generate force terms in child? Let's think. If local expansion is for acceleration directly, a constant term in parent acceleration polynomial? A constant acceleration (independent of position) has gradient? Actually local expansion for acceleration is a polynomial in η; β=0 term is constant acceleration. It should contribute to acceleration in child as constant, not zero. Our representation `a = Σ C_β ∇ η^β` cannot represent a constant acceleration term because gradient of η^β for β=0 is zero. We omitted constant acceleration. Does multipole far field produce constant acceleration term? For a single source, acceleration = -Gm/D^2 + 2Gmη/D^3 - ... has constant term. In our C_β representation, constant term came from β=1: C_1 ∇η^1 = C_1, so constant acceleration is represented by first-order coefficient, not β=0. Indeed ∇η^1 is constant. So β=0 term is always zero and not needed. A general polynomial acceleration can be represented by gradients of monomials; constant term from linear monomial. So L0 not used. In M2L formula, L0 corresponds to ∂^0? It would be unused. For translation L2L, linear terms etc. L0 not needed. But if we use potential coefficients, β=0 constant potential has no force but its derivatives for higher? No, potential constant no force. For acceleration coefficients, no β=0. Good. But M2L formula `L_beta = Σ_alpha M_alpha ∂^{α+β} S(D)/β!` for β=0 yields M0 S(D) (potential constant) not acceleration. Wait for β=1: M0 ∂S(D)/1! gives constant acceleration. So acceleration coefficient index β is order of local monomial whose gradient gives force. Thus β ranges 1..p. β=0 not needed. In M2L, for β from 1..p. Good. L2L translation for acceleration coefficients: Suppose parent local acceleration field around parent center: a(η_p) = Σ_{β} C_β ∇(η_p^β). Child center offset d, η_p = d + η_c. We need express as Σ C'_γ ∇(η_c^γ). This is just Taylor shift of the polynomial vector field. For pseudo code, `L2L(C_parent, d)` can be computed by evaluating polynomial and derivatives or using translation matrices. Since C_β are coefficients of a polynomial vector field, shift can be done with binomial/multinomial: If term C_β ∇ η_p^β. η_p^β = (d+η_c)^β = Σ_γ binomial(β,γ) d^{β-γ} η_c^γ. Gradient ∇_{η_p} η_p^β = ∇_{η_c} (d+η_c)^β? Since d constant, ∇_{η_p} η_p^β = ∇_{η_c} η_p^β = β_i (η_p)^{β-e_i}. Need express as gradient of η_c^γ times coefficients? There is identity: ∇ (d+η)^β = Σ_γ binomial(β,γ) d^{β-γ} ∇ η^γ? Let's test β=1: (d+η) gradient = 1 = ∇η^1, binomial β=1: γ=1 coefficient d^0, γ=0? binomial(1,0)=1 d^1 ∇η^0=0, so only γ=1. Good. β=2: ∇(d+η)^2 = 2(d+η) = 2d + 2η = 2d ∇η^1 + 2 ∇η^2? Since ∇η^1=1, ∇η^2=2η. Coeff: 2d from γ=1, 2 from γ=2. Formula Σ_γ multinomial(β,γ) d^{β-γ} ∇η^γ? For β=2, γ=1: binom(2,1)d^1 ∇η^1=2d; γ=2:1*∇η^2=2η. Good. β=0: ∇1=0. So translation: C'_γ += C_β * multinomial(β,γ) * d^{β-γ} for γ <= β (componentwise) maybe for each component? But C_β is scalar coefficient for gradient of scalar monomial. Since ∇η^β is vector. The identity holds vector-wise. So L2L: ``` for beta in 1..p: for gamma in 1..beta (multi-index <= beta): C_child[gamma] += C_parent[beta] * multinomial(beta, gamma) * d^(beta-gamma) ``` Where multi-index <= componentwise. This is for acceleration coefficients. Good. We can include this pseudo code. M2L formula with acceleration coefficients: Need ensure local coefficient C_β can be obtained by expanding source multipole acceleration field in η. We derived C_β = Σ_α M_α ∂^{α+β} S(D)/β! for β≥1? Let's test β=2: C_2 = M0 ∂^2 S/2! + M1 ∂^3 S/2!? For 1D source multipole M1=-mδ. Expansion of acceleration: G[-1/D^2 -2δ/D^3 -3δ^2/D^4 - ...]. Local terms: C1 ∇η + C2 ∇η^2 + C3 ∇η^3 ... C1 = M0 ∂S = -m/D^2 + M1 ∂^2S? Wait β=1: C1 = Σ_α M_α ∂^{α+1}S /1!. α=0: -m/D^2. α=1: M1 ∂^2S = -mδ * 2/D^3 = -2mδ/D^3. But expected constant acceleration -m/D^2 - 2δ/D^3. Good. So C1 includes constant and linear? C1 ∇η^1 = C1 constant. It includes both constant and? No C1 is constant; gradient of η^1 is 1. It includes -2δ/D^3 constant shift due to source offset. Good. C2 β=2: α=0 M0 ∂^2S/2! = (2/D^3)/2=1/D^3; α=1 M1 ∂^3S/2! = -mδ * (-6/D^4)/2 = 3mδ/D^4. L2P C2 ∇η^2 = 2C2 η = 2η/D^3 + 6δ η/D^4. Expected expansion of -1/(D+η+δ?)? Let's derive: source at δ, target η, distance D+η-δ? R = D + η - δ? Actually D=c_B-c_A, source δ, target η: r_i-r_j = D+η-δ. K=1/(D+η-δ). Acceleration ∂/∂η K = -1/(D+η-δ)^2. Expand in η: -1/D^2 (1 + η/D - δ/D)^-2 = -1/D^2 [1 + 2(η-δ)/D + 3(η-δ)^2/D^2 + ...] = -1/D^2 -2(η-δ)/D^3 -3(η-δ)^2/D^4 ... = -1/D^2 + 2δ/D^3 -2η/D^3 -3η^2/D^4 +6δη/D^4 ... Our local: C1 = -1/D^2 -2δ/D^3? Wait sign: M1 = (-1)^1 m δ = -mδ. ∂^2 S at D: 2/D^3. C1 includes M1*2/D^3 = -2δ/D^3. But expected constant +2δ/D^3? Let's check sign convention. Source at δ relative to A center, target at B center D away. Vector from source to target? R = r_i - c_A = D+η. Source position relative c_A = δ. Actual vector r_i - r_j = D+η-δ. K=1/|D+η-δ|. Taylor in δ: 1/(D+η-δ) = 1/(D+η) + δ/(D+η)^2 + ... (positive). Acceleration ∂η: -1/(D+η-δ)^2 = -1/(D+η)^2 - 2δ/(D+η)^3 + ... (negative δ). At η=0: -1/D^2 -2δ/D^3. So our C1 constant -1/D^2 -2δ/D^3 matches. Good. Expected sign was negative. Good. C2: from α=0 1/D^3, α=1 -mδ * ∂^3S/2. ∂^3 1/D in 1D = -6/D^4. product 3δ/D^4. C2=1/D^3+3δ/D^4. L2P 2C2η = 2η/D^3+6δη/D^4. Expected term linear in η: +2η/D^3? Expansion: -1/D^2 -2(η-δ)/D^3 -3(η-δ)^2/D^4... linear η: -2η/D^3? Wait compute carefully: Let u = η - δ. K = 1/(D+u) = 1/D - u/D^2 + u^2/D^3 - u^3/D^4 +... Acceleration = ∂η K = ∂u K = -1/D^2 + 2u/D^3 - 3u^2/D^4 + 4u^3/D^5... = -1/D^2 + 2(η-δ)/D^3 -3(η-δ)^2/D^4... Constant: -1/D^2 -2δ/D^3. Linear η: +2η/D^3 + 6δη/D^4? From -3(η^2 -2δη+δ^2)/D^4 linear: +6δη/D^4. Yes matches. Good. So formula works. Thus acceleration coefficient representation is coherent. Need M2P formula with M_α and D derivatives. We should define `M_α = Σ m (-1)^{|α|} δ^α / α!`. Then far acceleration for target in cell B from source cell A can be evaluated either M2P directly with R=r_i-c_A: `a = G Σ_α M_α ∇ D_α(R)`, where `D_α(R)=∂^α S(R)`. Good. M2L: For target cell B, from source A, local acceleration coefficients `C_B,β += G?` We can store without G and multiply at end. `C_β = Σ_{α+β<=p} M_α D_{α+β}(D) / β!`. Then L2P: `a = G Σ_{β>=1} C_β ∇ η^β`. M2L sum over all far source nodes A. L2L shift: `C_child,γ += Σ_{β>=γ} C_parent,β * multinomial(β,γ) * d^{β-γ}`. Need if C_parent includes G? We can store without G, apply G in L2P/P2P. Near P2P: direct acceleration with softening: ``` dx = xj - xi; dy = yj - yi r2 = dx*dx+dy*dy+eps2 inv_r3 = r2 ** -1.5 ax[i] += G * mj * dx * inv_r3 ``` This matches derivative of 1/r. Good. Integrator: ``` a0 = compute_accelerations() x += v*dt + 0.5*a0*dt^2 a1 = compute_accelerations() v += 0.5*(a0+a1)*dt ``` If storing mass, acceleration includes G*m_j; test mass cancels if acceleration. If force, multiply by m_i. Need mention periodic boundaries? If simulation in periodic box, kernel should be Ewald/periodic; FMM can be used with periodic image or Ewald. For simplicity, open boundary or non-periodic. If periodic, use periodic FMM or direct images; not in scope. Could mention. Need complexity: - Tree build O(N log N) or O(N) with sorting. - P2M O(N D K). - M2L O(N_far K^2) or O(N log N K^2) with well-separated lists. - L2P O(N K). - Near O(N Pmax K? Actually O(N Pmax)). With leaf capacity Pmax, O(N Pmax). Choose Pmax ~ 32-128. - Overall near-linear / O(N log N), memory O(N + Nnodes*K + farlists). Accuracy: - Truncation error per far interaction ~ O((h_s/h_t)^{p+1}); θ criterion: well-separated if `dist(c_A,c_B) > θ (h_A+h_B)` maybe; error controlled by θ^{p+1}. For θ=0.5, p=4 error ~0.03? Actually (0.5/1.5)^5 ~0.004. Good. Use p=6 for high accuracy. Softening avoids singular near; direct near handles close. Need maybe include "fast multiple method" could be "Fast Multipole Method" and "fast multiple" no? We'll say "I will assume Fast Multipole Method (FMM)." Now final answer structure: - Title. - Assumptions. - Approach: 1. Kernel and FMM math. 2. Spatial tree and well-separated criterion. 3. GPU mapping and kernels. 4. Integration and parameters. - Python pseudo code. - Notes for real implementation. Let's draft explanation in detail. Mathematical section: ``` Given N particles (r_i, m_i) in 2D. Acceleration: a_i = G Σ_{j≠i} m_j (r_j - r_i) / (|r_j-r_i|^2 + ε^2)^{3/2} ``` If true 2D log kernel: `a_i = 2πG Σ m_j (r_j-r_i)/(|...|^2+ε^2)`, replace S(r)=ln r. But we use 1/r. Quadtree node A: center c_A, half-width h_A. Multipole coefficients up to degree p: ``` M_A,α = Σ_{j∈A} m_j (-1)^{|α|} (r_j-c_A)^α / α! ``` Far evaluation M2P: ``` R = r_i - c_A a_i^far += G Σ_{|α|≤p} M_A,α ( ∂_x^{α_x+1}∂_y^{α_y} S(R), ∂_x^{α_x}∂_y^{α_y+1} S(R) ) ``` Local coefficients C_B,β from source A: ``` D = c_B - c_A C_B,β += Σ_{|α|+|β|≤p} M_A,α / β! * ∂^{α+β} S(D) ``` Then for particle in B: ``` η = r_i - c_B a_i^far += G Σ_{|β|≤p, β≠0} C_B,β ∇(η^β) ``` Near direct: ``` a_i^near += G Σ_{j in near leaves} m_j (r_j-r_i)/(r^2+ε^2)^{3/2} ``` M2L well-separated if boxes not in 3x3 neighborhood or `|c_A-c_B| > θ(h_A+h_B)`. GPU design: - SoA arrays on GPU: `x,y,vx,vy,m,ax,ay` float64/32. - Tree arrays: `cx,cy,h,parent,first,last` etc. `M,L` shape `(num_nodes, K)`. - Build/assign: compute leaf cell id from position, radix sort, count; if max occupancy > leaf_max, rebuild deeper. For adaptive quadtree, node assignment via occupancy. Rebuild every few steps. - Kernels: - `p2m`: parallel over particles/leaves; accumulate moments to ancestors. Use block reductions; atomic only final. - `m2l`: parallel over target nodes; each block loops far source list; uses shared memory for M; translation is small matrix multiply; can precompute translation tables for uniform grid. - `l2l`: parallel over nodes top-down; shift local coefficients to children. - `l2p`: one thread per particle; loads leaf local coefficients; evaluate K terms. - `near_p2p`: one block per leaf; load leaf+neighbor particles into shared memory; thread per particle; loop neighbors. - Use CUDA streams; no host sync; double precision if needed; mixed precision with correction; error control. Pseudo code: Need write in Python with comments. Let's make it coherent. Pseudo code skeleton: ``` # Python pseudo-code. Not a complete production code. import torch from torch import tensor from math import factorial from itertools import product G = 1.0 def multi_indices(p): return [(a,b) for a in range(p+1) for b in range(p+1-a)] def multi_factorial(alpha): return factorial(alpha[0]) * factorial(alpha[1]) def monomial(v, alpha): return v[0]**alpha[0] * v[1]**alpha[1] def grad_monomial(v, alpha): if alpha[0] > 0: gx = alpha[0]*v[0]**(alpha[0]-1)*v[1]**alpha[1] else: gx = 0.0 ... ``` `S_derivatives` placeholder: ``` def S_derivatives(R, max_order, eps): # returns 2D list D[a,b] # In real code: precomputed symbolic expressions for 1/sqrt(Rx^2+Ry^2+eps^2) # or evaluate recurrence using polynomial table. r2 = R[0]*R[0] + R[1]*R[1] + eps*eps D = [[0.0]*(max_order+1) for _ in range(max_order+1)] # pseudo: use automatic differentiation or unrolled formulas # D[a,b] = partial_x^a partial_y^b (1/sqrt(r2)) return D ``` But if we leave as placeholder, okay. Could include recurrence for Q in comments. Maybe better to include a pseudo function that computes via finite differences? No, not accurate. We can state "black box". Maybe include a small hardcoded for p=0? Not needed. QuadTree pseudo: ``` class QuadTree: def __init__(self, domain, max_depth, leaf_max): self.cx, self.cy, self.h, self.parent, self.children, self.first, self.last = ... def assign(self, x, y): # compute leaf id, sort, fill first/last def well_separated(self, a, b, theta): return abs(self.cx[a]-self.cx[b]) > theta*(self.h[a]+self.h[b]) and ... ``` Need box distance in 2D. Use center distance conservative. FMM class: ``` class NBody2D_FMM: def __init__(...): self.x = ... self.tree = QuadTree(...) self.p = p self.ALPHAS = multi_indices(p) self.K = len(self.ALPHAS) self.M = torch.zeros((tree.num_nodes, K), device='cuda', dtype=torch.float64) self.C = torch.zeros_like(self.M) # local acceleration coefficients self.far = tree.build_far_lists(theta) self.near = tree.build_near_leaf_lists() ``` P2M: ``` def p2m(self): self.M.zero_() # CUDA kernel pseudo (one thread per particle, atomic_add) for i in range(self.N): node = self.tree.leaf[i] while node >= 0: delta = (self.x[i]-self.tree.cx[node], self.y[i]-self.tree.cy[node]) for k, alpha in enumerate(self.ALPHAS): val = self.m[i] * ((-1)**(alpha[0]+alpha[1])) * monomial(delta, alpha) / multi_factorial(alpha) self.M[node,k].atomic_add(val) # pseudo; use block reduction in real code node = self.tree.parent[node] ``` But `torch.Tensor.atomic_add` not exists. Pseudocode okay. Could mark `# atomic_add in CUDA`. Maybe p2m bottom-up would be better in pseudo: ``` # leaf moments directly for leaf in tree.leaves: M[leaf] = sum over particles in leaf # internal nodes: use M2C translation from children for level reversed: for node: M[node] = 0 for child: M[node] += multipole_center_shift(M[child], center[child]-center[node], p) ``` This avoids atomics and is standard. Need provide `multipole_center_shift`? Could say black-box. But P2M direct to ancestors is simpler. However for FMM explanation, bottom-up with center shift is more proper. Could include both: pseudo uses direct ancestor accumulation for clarity; production uses bottom-up reductions. Let's do direct ancestor for pseudo but note production. M2L: ``` def m2l(self): self.C.zero_() # kernel: block per target node for b in range(num_nodes): for a in self.far[b]: D = (tree.cx[b]-tree.cx[a], tree.cy[b]-tree.cy[a]) Dmat = S_derivatives(D, self.p, self.eps) for kb, beta in enumerate(self.ALPHAS): if beta == (0,0): continue total = 0 for ka, alpha in enumerate(self.ALPHAS): if alpha[0]+alpha[1] + beta[0]+beta[1] <= self.p: total += self.M[a,ka] * Dmat[alpha[0]+beta[0], alpha[1]+beta[1]] / multi_factorial(beta) self.C[b,kb] += total # L2L top-down self.l2l() ``` But if `C` for all nodes includes far contributions, then L2L adds parent local to children; but if child already has its own far contributions, okay. Need not double count? Standard local expansion at child should include far fields from nodes far from child; parent local includes far fields from nodes far from parent. A source node far from parent may not be far from child? If parent and child boxes, if source well-separated from parent, it's also well-separated from child? Child is inside parent; distance from source to child maybe less, but if well-separated from parent using θ, likely well-separated from child? Maybe not if child near parent boundary and source adjacent to parent? But far list excludes near cells of parent; child near neighbor might include source that was near parent. Standard algorithm carefully avoids double counting by not M2L from nodes that are not well-separated. L2L propagation of parent L to child is for far fields that were well-separated from parent; child should not separately M2L from those same nodes if they are also well-separated from child? In standard FMM, M2L is done for well-separated cell pairs; then L2L propagates to descendants. A source cell may be well-separated from both parent and child; if M2L adds to both, double count? Standard algorithms use a specific list: for each cell, M2L from cells that are well-separated and not in the cell's near field, but L2L handles ancestor far fields; to avoid double counting, M2L for child should only include cells that are not well-separated from ancestor? Hmm. Need be careful. Standard FMM algorithm (quadtree): - P2M. - M2L: For each cell, sum M2L from all cells that are well-separated from it? I think yes, but then L2L? Let's recall. In tree FMM, local expansion at a cell is built from multipoles of well-separated cells. Then local expansions are passed to children (L2L). If child also receives M2L from same source, double. But standard M2L traversal for a target cell stops when a source node is well-separated from the target cell and adds M2L; it does not recurse into that source. If source is well-separated from parent, it may also be well-separated from child, so child's M2L would add it too if independent. However in standard algorithm, M2L is not done independently for every cell from scratch; it's done in a top-down manner? Let's recall from Greengard: For each cell, local expansion is computed from multipoles of cells that are not in the neighborhood. Then local expansion is translated to children. If you compute local expansion for every cell using all well-separated cells, you'd double count because ancestor local expansion includes same far cells? But maybe child's direct M2L only uses source cells that are well-separated from child but not well-separated from parent? The well-separated relation changes; a source far from parent may be not far from child? Actually if source is far from parent, distance > θ (h_parent + h_source). Child center is within h_parent/2 of parent center, child size h_parent/4. Minimum distance to child center could be distance - h_parent/2. Could still > θ(h_child+h_source)? For θ ~0.5, maybe yes. So double possible. Standard solution: M2L is done for pairs of cells that are well-separated and at the same "level" in a multi-level grid? In hierarchical, local expansion at a node is only from nodes that are well-separated and not in the ancestor path? Need precise. Simplify by avoiding L2L: For each particle, use M2P from all well-separated nodes (top-down) and direct near. This is unambiguous and correct, though not full FMM local. But user asked FMM; local expansion optional? FMM includes local expansions for efficiency, but a per-particle multipole evaluator is still using multipole expansions. However to be safe, we can present a simplified but correct FMM-like algorithm: use local expansions at leaf cells only, computed from all well-separated source nodes (no L2L). This avoids double counting if we define far lists per leaf and near lists per leaf. Complexity maybe O(N * far_nodes_per_leaf). If tree is adaptive, a leaf may have many far nodes but O(N log N). Correct: For each target leaf B, for every source node A not in near leaves of B and well-separated, add M2L to L_B. Then L2P for particles in B. No ancestor L. This is correct because we include all far source nodes. It may be less efficient than standard hierarchical L2L but still fast if far lists are not huge. For pseudo code, that's easier and avoids double count. We can mention production uses multi-level L2L to reduce far pair count. So pseudo: `m2l` computes local coefficients for each leaf (or node) from its `far_sources` list. `l2p` uses leaf C. No L2L. In explanation, mention standard full FMM uses M2L/L2L at all levels; simplified pseudo uses leaf-level far list for clarity. But user asked design; we can state production uses hierarchical local expansions; pseudo shows core. Need far_sources for leaf: all nodes that are well-separated from the leaf's 3x3 neighborhood? If a source node overlaps any of the 9 neighbor leaves, it's near and handled by direct P2P. If not, and well-separated, use multipole. If a source node is ancestor of leaf? It contains target and near; skip. If source node is descendant of leaf? Leaf has no descendants; if target leaf, source leaf same/near direct; other source leaves not descendants. For internal nodes, far lists maybe for leaves only. Good. Near leaves: For each leaf, 8 neighboring leaves plus itself at same tree depth? In adaptive tree, neighboring cells may have different depths. Standard near field: for a particle, direct interact with all particles in the 3x3 neighborhood of its leaf at the leaf level. If neighbor cell is subdivided (finer), include all its subcells? If using leaf cells (maximal empty? all leaves same depth if fixed tree). Simplify: fixed-depth quadtree (uniform grid of cells) with leaf capacity. Then near leaves are 9 grid cells. This is easiest. But fixed-depth uniform grid with local M2L from all far cells is O(N^2) if many cells? Let's analyze: C^2 cells, each far list O(C^2) => M2L O(C^4) = O(N^2) if N~C^2. Not good. But if C ~ N^{1/4}? To keep leaf capacity constant, C ~ sqrt(N) (since 2D cells C^2, N/C^2 = P). Then far pairs per cell ~C^2, total C^4 ~ N^2. Not fast. Multi-level needed. So fixed-depth not acceptable for large N. Need hierarchical. Could use quadtree with leaf-level far lists; number far pairs O(N log N)? Let's check. For each leaf, far nodes include many larger cells. Number of well-separated pairs between leaves and all nodes maybe O(N log N). Could be okay. If we compute M2L for each leaf from all far nodes, total far pairs maybe O(N log N) or O(N^2)? In a quadtree with N leaves, number of well-separated pairs between leaves and nodes (including internal) might be O(N log N). Standard FMM uses local expansions at all levels to get O(N). But O(N log N) is acceptable? Maybe. For pseudo, okay. Alternatively, present standard full FMM with M2L and L2L but not get bogged in double count; we can state "the far list is built so each well-separated source/target pair is counted once in the hierarchical scheme (standard FMM bookkeeping)". That's acceptable for high-level design. Need maybe avoid overly specific incorrect algorithm. Let's formulate explanation at design level, not exact low-level bookkeeping. Pseudo code can be simplified with comments "standard FMM traversal". For correctness, we can define far list per target node from a standard well-separated traversal that does not double count with L2L. But in pseudo code we can just use `far_sources[b]` precomputed by a correct FMM list builder. That's fine. Let's think about standard FMM bookkeeping to not state wrong. In hierarchical FMM, one algorithm: - For each cell C, local expansion L_C is computed from multipoles of cells that are well-separated from C and not in the same "neighbor" list. This can be done by a top-down traversal: function `M2L(target)` recursively visits source tree: if source node is well-separated from target, do M2L(source,target); else if source is ancestor of target? skip; else recurse into children. If target is a leaf, this will add M2L from all well-separated source nodes. Then L2L from parent to child? If we do this for every leaf, double count as noted. But if we do M2L for every node (not just leaves) and then L2L, we must ensure M2L for a node only from source nodes that are well-separated from node but not from its ancestors? Actually standard: local expansion at a node is built from multipoles of well-separated nodes at the same or lower level? Let's search memory: In FMM quadtree, algorithm steps: P2M, M2L, L2P, L2L. M2L: For each cell, for each well-separated cell, translate multipole to local. L2L: top-down translate local from parent to children. If M2L is done for all cells from all well-separated cells, L2L will double count. But perhaps M2L only considers cells that are well-separated and *not* well-separated from the parent? Or M2L is done before L2L and local expansion of a cell includes contributions from cells that are not in its near field, including those far from parent; then L2L adds contributions that are far from parent but not necessarily far from child? Wait if parent local already includes far from parent, child should inherit those. Child direct M2L should not include those same far-from-parent cells. How to know? The M2L traversal for child can stop when a source node is well-separated from the *target's leaf*, but if it was already included in ancestor local, it would double. The standard way is to perform M2L in a top-down manner where local expansion from parent is propagated, and M2L for a cell only adds source nodes that are well-separated from the cell but were not already accounted for by ancestors. This is often done by using a list of "well-separated" nodes that are not in the ancestor's far list? Complex. Maybe use multi-level uniform grid FMM where near field interactions handled at finest level, and far field at each level with well-separated cells; no L2L? There is an algorithm (M3) where at each level, compute M2L for cells at that level from cells well-separated at that level, and L2P; plus near P2P at finest. Multi-level avoids double count by using different levels? Could be complex. Given this is a design answer, not production code, we can abstract the far field stage: "Use a standard tree traversal to build local expansions; the exact bookkeeping is the same as Barnes-Hut/FMM: for each particle, a far field node is used if it is well separated, otherwise descend; local expansions are cached per cell." Then pseudo code can be conceptual. The user likely doesn't require exact FMM list algorithm. Maybe a simpler correct GPU approach: per-particle tree traversal with multipole M2P and near P2P. This is basically Barnes-Hut, not full FMM but uses multipole expansion. It is easy to make correct: For each target particle i, traverse tree: ``` def force(i): a=0 stack=[root] while stack: node=pop if node == leaf(i) or not well_separated(node, leaf_i): if node is leaf: direct P2P with particles in node else: push children else: a += M2P(node, r_i) return a ``` This counts each target's force once. It uses multipole for far nodes. It doesn't use local expansion, but it's a "fast multipole" in loose sense? Fast Multipole Method specifically uses local expansions to reduce per-particle far cost. But per-particle traversal with M2P is O(N log N) and GPU-friendly? It may be okay for moderate N. But user explicitly "Use the fast multiple method" likely expects FMM stages. We can include local expansion as an optimization but pseudo code can be per-particle M2P for clarity? Need decide. Maybe present both: "Core design: hierarchical FMM with M2L/L2P. For readability, the pseudo-code shows a simplified per-particle multipole evaluator; replace `evaluate_far` with cached local expansions in production." But the user asked return explanation and pseudo code; a simplified pseudo code is acceptable if explanation notes full FMM. Let's aim for a balanced answer: Explanation describes full FMM (P2M, M2L, L2P, near). Pseudo code implements a conceptually correct simplified version using leaf-level local expansions and a far list precomputed by standard well-separated traversal. It won't delve into double-count bookkeeping but states `far_sources` is the standard FMM far list (each source used once). That should satisfy. Need maybe include "well-separated" definition for quadtree: For two nodes A, B, if their bounding boxes do not touch and the minimum distance between boxes > θ * (size_A + size_B), then well-separated. More conservative: use center distance > θ*(h_A+h_B). In tree traversal, a node is far if it is not in the 3x3 neighborhood of the target cell. For variable depth, use box distance. Pseudo code details: - Use `torch.cuda` tensors. But pseudo loops in Python are not GPU; mark "CUDA kernel would replace loops with `for` over threads". Use comments like `# CUDA: one thread per i`. - Could include `@cupy`? Maybe PyTorch more common. Let's write pseudo code in a clear way: ``` # -------------------------------------------------------------- # Python pseudo-code for a GPU 2D N-body FMM # Not production code: loops are conceptual GPU work-items. # -------------------------------------------------------------- import torch from math import factorial G = 1.0 def multi_indices(p): return [(a,b) for a in range(p+1) for b in range(p+1-a)] def mf(alpha): return factorial(alpha[0]) * factorial(alpha[1]) def monomial(v, alpha): return v[0] ** alpha[0] * v[1] ** alpha[1] def grad_monomial(v, alpha): ax = alpha[0] * v[0] ** (alpha[0]-1) * v[1] ** alpha[1] if alpha[0] else 0.0 ay = alpha[1] * v[0] ** alpha[0] * v[1] ** (alpha[1]-1) if alpha[1] else 0.0 return ax, ay def inv_r_derivatives(R, max_order, eps): """ Return D[a,b] = partial_x^a partial_y^b (1/sqrt(Rx^2+Ry^2+eps^2)) for 0 <= a+b <= max_order. In real CUDA: precompute symbolic polynomials for p=4..8 or use AD. """ # pseudo black-box ... ``` Need maybe implement placeholder with zeros? It's pseudo. Could include a note: "For p<=4, unroll formulas; for p<=8, generate offline with sympy." QuadTree pseudo: ``` class QuadTree: def __init__(self, xmin,ymin,xmax,ymax, max_depth, leaf_max): self.cx = torch.zeros(max_nodes, device='cuda', dtype=torch.float64) ... self.parent = torch.full(...) self.child = torch.full((max_nodes,4), -1, dtype=torch.int32) self.first = torch.full((max_nodes), 0, dtype=torch.int32) self.last = torch.full((max_nodes), 0, dtype=torch.int32) self.leaf = torch.zeros(N, dtype=torch.int32) def assign(self, x, y): # CUDA: compute leaf cell for each particle, sort by leaf, count, # split until occupancy <= leaf_max, build node arrays. ... def well_separated(self, a, b, theta): # box distance conservative dx = abs(self.cx[a]-self.cx[b]) - (self.h[a]+self.h[b]) dy = abs(self.cy[a]-self.cy[b]) - (self.h[a]+self.h[b]) if dx < 0 or dy < 0: # boxes overlap/touch? maybe need 3D box distance return False return (dx*dx + dy*dy) ** 0.5 > theta * (self.h[a]+self.h[b]) ``` Box distance in 2D: If intervals overlap, distance 0. Conservative center distance is simpler: ``` return hypot(cx[a]-cx[b], cy[a]-cy[b]) > theta * (h[a]+h[b]) ``` Use that. Build far lists: ``` def build_far_lists(self, leaves, theta): far_ptr = torch.zeros(len(leaves)+1, dtype=torch.int32) # standard FMM traversal: for each target leaf, collect source nodes # that are well-separated and not in the target's 3x3 near neighborhood. # Store CSR. In production, build once per tree on GPU using tree traversal. ``` Need not implement. FMM class pseudo: ``` class NBodyFMM: def __init__(self, x, y, m, max_depth=8, p=4, theta=0.7, eps=1e-3, leaf_max=64): self.x = torch.tensor(x, dtype=torch.float64, device='cuda') ... self.tree = QuadTree(...) self.tree.assign(self.x, self.y) self.ALPHAS = multi_indices(p) self.K = len(self.ALPHAS) self.M = torch.zeros((tree.num_nodes, K), dtype=torch.float64, device='cuda') self.L = torch.zeros_like(self.M) self.ax = torch.zeros_like(self.x) self.ay = torch.zeros_like(self.y) self.far = self.tree.build_far_lists(theta) self.near = self.tree.build_near_leaf_lists() ``` P2M: ``` def p2m(self): self.M.zero_() # CUDA kernel: one thread per particle; depth small. for i in range(self.N): node = self.tree.leaf[i] while node >= 0: delta = (self.x[i]-self.tree.cx[node], self.y[i]-self.tree.cy[node]) for k, alpha in enumerate(self.ALPHAS): # atomic_add in real CUDA self.M[node,k] += self.m[i] * ((-1)**(alpha[0]+alpha[1])) * monomial(delta, alpha) / mf(alpha) node = self.tree.parent[node] ``` Could note this is simple; production uses leaf block reductions and bottom-up M2C. M2L: ``` def m2l(self): self.L.zero_() # CUDA kernel: one block per target leaf/node for b in self.tree.leaf_ids: # target cells needing local expansion for a in self.far[b]: D = (self.tree.cx[b]-self.tree.cx[a], self.tree.cy[b]-self.tree.cy[a]) Dtab = inv_r_derivatives(D, self.p, self.eps) for kb, beta in enumerate(self.ALPHAS): if beta == (0,0): continue s = 0.0 for ka, alpha in enumerate(self.ALPHAS): if alpha[0]+alpha[1] + beta[0]+beta[1] <= self.p: s += self.M[a,ka] * Dtab[alpha[0]+beta[0], alpha[1]+beta[1]] / mf(beta) self.L[b,kb] += s # optional standard L2L top-down if far lists are hierarchical # self.l2l() ``` Maybe include L2L optional: ``` def l2l(self): # top-down shift local polynomial from parent to children for level in range(self.max_depth): for node in nodes_at_level(level): for child in self.tree.child[node]: if child == -1: continue d = (self.tree.cx[child]-self.tree.cx[node], ...) newL = torch.zeros(K) for kb,beta in enumerate(ALPHAS): for kg,gamma in enumerate(ALPHAS): if gamma <= beta componentwise: coeff = multinomial(beta,gamma) * monomial(d, beta-gamma) newL[gamma] += L[node,kb]*coeff L[child] += newL ``` But if we call L2L after M2L for all nodes, double count concerns. Could mark as optional production bookkeeping. Maybe not include in pseudo to keep simple. L2P: ``` def l2p(self): # one thread per particle for i in range(self.N): b = self.tree.leaf[i] eta = (self.x[i]-self.tree.cx[b], self.y[i]-self.tree.cy[b]) for k, beta in enumerate(self.ALPHAS): if beta == (0,0): continue gx, gy = grad_monomial(eta, beta) self.ax[i] += G * self.L[b,k] * gx self.ay[i] += G * self.L[b,k] * gy ``` Near P2P: ``` def near_p2p(self): # one block per leaf; load neighbor particles in shared memory for b in self.tree.leaf_ids: nb_list = self.near[b] # 9 leaves, or all neighboring leaf cells # load particles in union nb_list into shared memory for i in self.tree.range(b): totalx=totaly=0 for nb in nb_list: for j in self.tree.range(nb): if i == j: continue dx = self.x[j]-self.x[i] dy = self.y[j]-self.y[i] r2 = dx*dx + dy*dy + self.eps2 inv_r3 = r2 ** -1.5 totalx += self.m[j] * dx * inv_r3 totaly += self.m[j] * dy * inv_r3 self.ax[i] += G * totalx self.ay[i] += G * totaly ``` Compute accelerations: ``` def compute_accelerations(self): self.ax.zero_(); self.ay.zero_() self.p2m() self.m2l() self.l2p() self.near_p2p() ``` Step: ``` def step(self, dt): a0x,a0y = self.ax.clone(), self.ay.clone() self.x += self.vx*dt + 0.5*a0x*dt*dt self.y += self.vy*dt + 0.5*a0y*dt*dt self.compute_accelerations() self.vx += 0.5*(a0x + self.ax)*dt self.vy += 0.5*(a0y + self.ay)*dt ``` Need if tree static, positions moved; need reassign/rebuild after step or every k. Add: ``` if step % rebuild_every == 0: self.tree.assign(...) ``` But velocity Verlet needs a at new positions; if tree rebuilt between a0 and a1? Usually compute a1 with updated positions and same tree? Tree should be updated before a1 or periodic. Pseudo: ``` def step(dt): a_old = compute_accelerations() update positions rebuild_tree_if_needed() a_new = compute_accelerations() update velocities ``` Good. Need include `multinomial` if L2L. Could skip. Maybe include "If you want true 2D logarithmic kernel" in pseudo? Could add note: replace `S(R)=log(r)` and derivatives; moments similar. But not in code. Potential improvements: - For M2L, translation matrices can be precomputed for relative node center offsets in uniform grids; in quadtree, use formulas. Use `torch.matmul` per block? In CUDA, use shared memory and unrolled loops. - Use double precision for moments/local to avoid cancellation; near direct can use double. If using single precision, use double for tree moments or error correction. - Use softening `eps` and skip self. - Use `leaf_max` to control near cost and accuracy; if a cell overflows, subdivide. - Use `theta` and `p` tradeoff: smaller θ or larger p more accurate, slower. Need maybe include "NVIDIA GPU" specifics: Use CUDA, shared memory, coalesced access, `cudaMemcpyAsync`, streams, `cupy`/PyTorch. Mention `nvidia gpu` not specific. Could mention `CuPy` pseudo code maybe more GPU-like. PyTorch loops pseudo okay. Maybe include a note: "For a real NVIDIA implementation, the kernels marked 'CUDA kernel' would be written in CUDA C++ or generated by Numba CuPy; the Python pseudo-code expresses the parallel work distribution." Good. Let's think about possible user expectation: "fast multiple method" maybe "Fast Multiple Method"? Could be "fast multipole method". Need maybe not overemphasize typo. Say "I assume Fast Multipole Method (FMM)." Good. Need ensure final answer not too terse. Desired oververbosity 9, so detailed. But not insane. We can include equations and pseudo code. Let's refine mathematical formulas for local expansion with acceleration coefficients. Define multi-index α = (α_x, α_y), |α|=α_x+α_y. `α! = α_x! α_y!`, `R^α = R_x^{α_x}R_y^{α_y}`. Kernel `K(r)=1/√(r_x^2+r_y^2+ε^2)`. Acceleration: `a_i = G Σ_j m_j ∇_{r_i} K(r_i-r_j)`. Multipole coefficients: `M_A,α = Σ_{j∈A} m_j (-1)^{|α|} (r_j-c_A)^α / α!`. M2P: `a_i(A) = G Σ_{|α|≤p} M_A,α ∇ ∂^α K(R_i)` where `R_i = r_i - c_A`, and `∂^α = ∂_x^{α_x}∂_y^{α_y}`. Local coefficients for cell B from A: `C_B,β += Σ_{|α|+|β|≤p} M_A,α / β! * ∂^{α+β} K(c_B-c_A)`. L2P: `a_i(B) = G Σ_{|β|≤p, β≠0} C_B,β ∇(r_i-c_B)^β`. Check dimensions: C_B,β has units mass/length^{?}. For β=1, C1 mass/length^2; ∇η^1 dimensionless -> acceleration mass/length^2 times G. Good. M2L formula uses `∂^{α+β} K(D)`. If β=(1,0), α=0 -> ∂xK. Good. M2P and M2L both use `K` derivatives. Need derivative up to p for M2L, p+1 for M2P if evaluating directly. In local approach, L2P uses polynomial derivatives not K derivatives. M2L uses up to p. If we also use M2P fallback, p+1. Near direct with softening: `a_i += G m_j (r_j-r_i) / (|r_j-r_i|^2+ε^2)^{3/2}`. This is consistent with ∇K. Well-separated criterion: Use boxes `[c_x-h_x, c_x+h_x] x [c_y-h_y, c_y+h_y]`. Let `d_box` be minimum distance between boxes. If `d_box > θ (h_A+h_B)` or more simply `|c_A-c_B| > θ(h_A+h_B)`, far. Need choose θ ~ 0.5-1.0. For accuracy, use `d_box` not center distance? Center distance can be conservative (more far? Actually if boxes far, center distance > box distance; using center distance > θ(h_A+h_B) may classify some overlapping? If boxes overlap, center distance may still > θ sum? Example same size h, adjacent boxes center distance 2h, θ=0.5 => 2h > h true but boxes touch? Need exclude neighbors. So need box distance or 3x3 neighborhood. Use box distance: if boxes intersect or touch, not far. If `d_box > θ (h_A + h_B)`. For adjacent boxes d_box=0, not far. For diagonal adjacent d_box=0? If touching at corner, 0. Good. For farther, distance positive. Use this. In pseudo `well_separated`: ``` dx = max(0, abs(cx[a]-cx[b]) - h[a] - h[b]) dy = max(0, abs(cy[a]-cy[b]) - h[a] - h[b]) return math.hypot(dx, dy) > theta * (h[a]+h[b]) ``` This is minimum box distance. Good. Near leaves: Direct for nodes whose boxes are not well-separated. For a target leaf, near source leaves are those in its 3x3 neighborhood (including diagonal) at the same depth. For adaptive, include any source leaf whose box is not well-separated from target leaf. In pseudo, `near[b]` built by checking all leaves? O(C^2) but once. Could say standard neighborhood list. Need maybe include "double counting" in near: If `near[b]` includes both b and nb, and we compute acceleration for each i against all j in near leaves, for pair i in b, j in nb, i's acceleration includes j; j's acceleration computed in j's block includes i. That's correct. If using Newton's third law to halve work, need symmetry; not in pseudo. Need maybe include "mass units" acceleration vs force. We'll compute acceleration (divide by m_i omitted). If user wants force, multiply by m_i. Potential issue: In M2L, if local coefficients C_B,β are summed over far source nodes A, and near direct handles near nodes, total acceleration = far + near. Good. Need mention periodic boundaries: If periodic, direct near should include images and FMM should be periodic or use Ewald; not covered. Could add at end. Let's consider if 2D N-body gravity should use logarithmic potential. In 2D, Poisson equation ∇²Φ=2πGρ gives Φ=(G/2π?) log r? Actually force = ∇²? In 2D, gravitational force between line masses is proportional 1/r, potential log r. Many "2D N-body" simulations still use 3D 1/r with x,y only. Need mention: "If by 2D you mean a strictly 2D Poisson/log kernel, replace K(r)=1/r by K(r)=ln r (or -ln r depending sign) and use the same FMM structure; the pseudo-code's `inv_r_derivatives` becomes log-derivatives." Good. Need maybe include "fast multiple method" could be "fast multipole" and "multiple" maybe "multipole". Good. Let's think about pseudo code's `inv_r_derivatives`. Could we provide a concrete recurrence for D values to make code more complete? Let's try to derive a numerical recurrence for D[a,b] without polynomial derivatives maybe using previous D values? There is a recurrence for derivatives of r^{-1}: `∂_x^{a+1}∂_y^b f = -x ∂_x^a ∂_y^b f - r^2 ∂_x^{a+2} ∂_y^b f`? Let's test. Use Euler: ∂_x (x f) = f + x ∂_x f = f - x^2 r^{-3}. Not. Maybe recurrence from `∂_x f = -x f / r^2`. Then for D[a,b] = ∂_x^a∂_y^b f. For a>=1: D[a,b] = ∂_x^{a-1}∂_y^b (∂_x f) = ∂_x^{a-1}∂_y^b (-x f/r^2). Hard. Could use relation: `∂_x^{a+1} f = -x ∂_x^a f - y^2 ∂_x^a f - ε^2 ∂_x^a f + ∂_x^a(∂_x(???))` Not. Let's derive from Q recurrence but need Q derivatives. We can compute Q polynomials symbolically for p small. In pseudo, okay. Maybe include a simple pseudo implementation using finite differences? No. Could use complex analytic functions? For 1/r in 2D, derivatives can be computed via recurrence in polar coordinates. But Cartesian local expansion easier with black-box. Let's include a note: "The only numerically delicate part is generating derivatives of `1/r`. For p=4 or 6, generate them once with sympy or hard-code. The recurrence for the polynomial Q is: Q[0,0]=1, Q[a,b]=r² ∂_x Q[a-1,b] - (2(a+b)-1)x Q[a-1,b], etc." This gives implementer a path. Let's write it in explanation or pseudo comment. Detailed recurrence: Let `D_{a,b} = r^{-(2(a+b)+1)} Q_{a,b}`. Base `Q_{0,0}=1`. For `a>0`: `Q_{a,b} = r^2 ∂_x Q_{a-1,b} - (2(a+b)-1) x Q_{a-1,b}`. For `b>0`: `Q_{a,b} = r^2 ∂_y Q_{a,b-1} - (2(a+b)-1) y Q_{a,b-1}`. Need `r^2 = x^2+y^2+ε^2`. This computes Q polynomial. To evaluate D, need Q and maybe Q derivatives for next. We can store Q polynomials symbolically. Good. For gradient of ∂^α K, just use D[α+e]. So need D up to p+1. Could include in pseudo: ``` # In production, generate Q[a,b] polynomials up to p+1 offline. # Then D[a,b] = Q[a,b](Rx,Ry) / r2**(a+b+0.5) ``` Good. Now, think about GPU performance for M2L translation. The translation matrix size K x K? K=(p+1)(p+2)/2. For p=4 K=15, matrix 225; p=6 K=28, 784. Small. In CUDA, one warp could handle one source-target pair, threads over β/α. Use shared memory. Precompute `partial_S(D, α+β)` for each pair? It depends on D. For uniform grid, D takes finite values; precompute translation tables per relative offset. For quadtree, compute on the fly. Could mention. Tree building details: - Fixed domain with maximum coordinate. If particles leave domain, re-root or use periodic. For open boundaries, domain must contain particles; rebuild with new bbox. - Adaptive quadtree: start root, subdivide if count > leaf_max or depth < max_depth. To do on GPU: recursive subdivision not easy; use grid-based at max depth then compress? For pseudo, assume CPU builds tree or GPU uses max-depth grid. In production, rebuild on GPU by counting max-depth cells, then building a tree of occupied cells using prefix sums. Or rebuild every N steps on CPU and copy. Mention. - For dynamic simulation, moving particles may change tree; FMM can use a static tree if particles don't move far relative to cell size; rebuild when needed. Or use "moving particle FMM" with local updates. For pseudo, rebuild each step or every few steps. Need maybe include memory estimate: - Particles: 6 doubles ~48 bytes each. - Nodes: moments 2*K doubles. For p=4 K=15 -> 240 bytes/node. If nodes ~N, ~240N bytes plus tree ~ maybe. Far lists. Good for GPU memory. Potential edge cases: - Self interaction: skip `i==j` in near; M2P/M2L doesn't include target node if not well-separated; ensure far lists exclude ancestors and near cells. - Softening: Use same ε in derivatives? `K(r)=1/√(r²+ε²)`. Derivatives formulas include ε. - Overflows/underflows: Use double; if r2 small in far? Far distance large, okay. Near direct with softening. - Accuracy validation: Compare to direct O(N²) for small N, check energy drift, convergence with p, θ, leaf_max. Now, let's decide final wording. Start: "Below I assume 'fast multiple method' means Fast Multipole Method (FMM). I design a 2D quadtree FMM for Newtonian 1/r gravity in the x-y plane. If you truly need 2D log gravity, swap kernel." Then sections. Let's write explanation in final: ### 1. Physical model Equation. Softening. Acceleration. ### 2. FMM representation Define multi-index, kernel, multipole, local. Provide equations. ### 3. Tree and well separated Quadtree, leaf capacity, near 3x3, far criterion. ### 4. GPU mapping SoA, kernels, parallelization, shared memory, atomics, rebuild, precision. ### 5. Python pseudo code Large code block. ### 6. Practical parameters and validation p, θ, leaf_max, dt, validation. Need ensure pseudo code not too long but enough. Let's craft code carefully. Pseudo code code block: ``` # Python pseudo-code. The loops are not Python execution loops; # they indicate CUDA thread/block work. Use CuPy/PyTorch/CUDA for real code. import torch from math import factorial, hypot G = 1.0 def multi_indices(p): return [(a,b) for a in range(p+1) for b in range(p+1-a)] def mf(alpha): return factorial(alpha[0]) * factorial(alpha[1]) def monomial(v, alpha): return v[0] ** alpha[0] * v[1] ** alpha[1] def grad_monomial(v, alpha): ax = alpha[0] * v[0] ** (alpha[0]-1) * v[1] ** alpha[1] if alpha[0] else 0.0 ay = alpha[1] * v[0] ** alpha[0] * v[1] ** (alpha[1]-1) if alpha[1] else 0.0 return ax, ay def inv_r_derivatives(R, max_order, eps): # Black box: D[a,b] for 0<=a+b<=max_order. # Generate Q polynomials offline: # D[a,b] = r**(-(2*(a+b)+1)) * Q[a,b](Rx,Ry) # Q[0,0]=1 # Q[a,b] = r2*dQdx(a-1,b) - (2*(a+b)-1)*Rx*Q[a-1,b] # Q[a,b] = r2*dQdy(a,b-1) - (2*(a+b)-1)*Ry*Q[a,b-1] # r2 = Rx*Rx + Ry*Ry + eps*eps # For p=4/6, unroll the formulas. D = [[0.0]*(max_order+1) for _ in range(max_order+1)] # pseudo: fill D return D ``` Maybe the list dimensions max_order+1 for each a,b but only a+b<=max_order used. QuadTree: ``` class QuadTree: def __init__(self, domain, max_depth, leaf_max): xmin,ymin,xmax,ymax = domain self.max_depth = max_depth self.leaf_max = leaf_max # node arrays on GPU (pseudo) self.cx = torch.zeros(1, dtype=torch.float64, device='cuda') ... def assign(self, x, y): # 1) compute candidate max-depth cell for each particle # 2) radix sort by cell # 3) build quadtree nodes from occupied cells, subdivide until count <= leaf_max # 4) store leaf[i], parent, child, first, last pass def well_separated(self, a, b, theta): dx = max(0.0, abs(self.cx[a]-self.cx[b]) - self.h[a] - self.h[b]) dy = max(0.0, abs(self.cy[a]-self.cy[b]) - self.h[a] - self.h[b]) return hypot(dx, dy) > theta * (self.h[a] + self.h[b]) def build_near_lists(self): # for each leaf, leaves in its 3x3 neighborhood (or all non-far leaves) pass def build_far_lists(self, theta): # standard FMM traversal; CSR far_sources per target leaf # exclude ancestors, self, and near neighbors pass ``` Need maybe `num_nodes` property. FMM class: ``` class NBody2D_FMM: def __init__(self, x, y, m, domain, max_depth=8, p=4, theta=0.7, eps=1e-3, leaf_max=64): self.x = torch.tensor(x, dtype=torch.float64, device='cuda') ... self.tree = QuadTree(domain, max_depth, leaf_max) self.tree.assign(self.x, self.y) self.ALPHAS = multi_indices(p) self.K = len(self.ALPHAS) self.M = torch.zeros((self.tree.num_nodes, self.K), dtype=torch.float64, device='cuda') self.L = torch.zeros_like(self.M) self.ax = torch.zeros_like(self.x) self.ay = torch.zeros_like(self.y) self.far = self.tree.build_far_lists(theta) self.near = self.tree.build_near_lists() ``` P2M with atomic add pseudo: ``` def p2m(self): self.M.zero_() # CUDA: one thread per particle; use atomicAdd(double) or block reduction. for i in range(self.N): node = self.tree.leaf[i] while node >= 0: dx = self.x[i] - self.tree.cx[node] dy = self.y[i] - self.tree.cy[node] for k, (a,b) in enumerate(self.ALPHAS): sign = -1.0 if (a+b) % 2 else 1.0 val = sign * self.m[i] * (dx**a) * (dy**b) / mf((a,b)) # atomic_add(self.M[node,k], val) self.M[node,k] = self.M[node,k] + val node = self.tree.parent[node] ``` Potential issue: internal node moments about node center computed by summing particle contributions to each ancestor. Correct. But if a particle is in leaf, its ancestors include root; good. Complexity O(N depth K). M2L: ``` def m2l(self): self.L.zero_() # CUDA: one block per target leaf. far is CSR. for b in self.tree.leaf_ids: for a in self.far[b]: D = (self.tree.cx[b]-self.tree.cx[a], self.tree.cy[b]-self.tree.cy[a]) Dtab = inv_r_derivatives(D, self.p, self.eps) for kb, beta in enumerate(self.ALPHAS): if beta == (0,0): continue s = 0.0 for ka, alpha in enumerate(self.ALPHAS): if alpha[0]+alpha[1] + beta[0]+beta[1] <= self.p: s += self.M[a,ka] * Dtab[alpha[0]+beta[0], alpha[1]+beta[1]] / mf(beta) self.L[b,kb] = self.L[b,kb] + s ``` Could include `if not tree.well_separated(a,b): continue` but far list prechecked. Near P2P: ``` def near_p2p(self): # CUDA: one block per leaf; load neighbor particles to shared memory. for b in self.tree.leaf_ids: nb_list = self.near[b] # optional: build local arrays sx,sy,sm for i in range(self.tree.first[b], self.tree.last[b]): ax = 0.0; ay=0.0 for nb in nb_list: for j in range(self.tree.first[nb], self.tree.last[nb]): if i == j: continue dx = self.x[j] - self.x[i] dy = self.y[j] - self.y[i] r2 = dx*dx + dy*dy + self.eps2 inv_r3 = r2 ** -1.5 ax += self.m[j] * dx * inv_r3 ay += self.m[j] * dy * inv_r3 self.ax[i] = self.ax[i] + G * ax self.ay[i] = self.ay[i] + G * ay ``` Need `range(first,last)` if first/last inclusive? Use pseudo. L2P: ``` def l2p(self): for i in range(self.N): b = self.tree.leaf[i] eta = (self.x[i]-self.tree.cx[b], self.y[i]-self.tree.cy[b]) ax=0.0; ay=0.0 for k,beta in enumerate(self.ALPHAS): if beta == (0,0): continue gx, gy = grad_monomial(eta, beta) ax += self.L[b,k] * gx ay += self.L[b,k] * gy self.ax[i] = self.ax[i] + G*ax self.ay[i] = self.ay[i] + G*ay ``` Order: `l2p` then `near_p2p` or vice versa; both add. Compute: ``` def compute_accelerations(self): self.ax.zero_(); self.ay.zero_() self.p2m() self.m2l() self.l2p() self.near_p2p() ``` Step: ``` def step(self, dt, rebuild_every=10): a0x,a0y = self.ax.clone(), self.ay.clone() self.x += self.vx*dt + 0.5*a0x*dt*dt self.y += self.vy*dt + 0.5*a0y*dt*dt if self.step_count % rebuild_every == 0: self.tree.assign(self.x,self.y) self.far = self.tree.build_far_lists(self.theta) self.near = self.tree.build_near_lists() self.compute_accelerations() self.vx += 0.5*(a0x + self.ax)*dt self.vy += 0.5*(a0y + self.ay)*dt self.step_count += 1 ``` Need initialize `step_count=0`. Maybe include `compute_accelerations` before first step. Add note: In real code, `M` and `L` could be allocated once max size; tree node count may change; use dynamic allocation or max_nodes. Potential bug: In P2M, using `self.M[node,k] = self.M[node,k] + val` in a GPU kernel is race; pseudo okay with atomic. Use comment. Potential bug: M2L uses local coefficients for target leaf only. If far list standard hierarchical includes only some sources and L2L not used, may miss. But pseudo says far list is standard FMM far list that covers all far contributions. Good. Maybe include optional `m2p_fallback` for clarity? Not necessary. Need mention "The local expansion coefficients `L` are for acceleration, not potential. If you store potential coefficients, take gradient with sign." Good. Let's consider if `M2L` formula with `Dtab[alpha+beta] / mf(beta)` is enough for acceleration coefficients. We derived yes. But need `β!` denominator. In code `mf(beta)`. Good. Need maybe include `alpha` and `beta` multi-index ordering. Good. Need maybe mention that for `p` degree, local expansion should include β up to p; but in M2L sum α+β <= p. If we want more accurate, use p_M and p_L maybe. Good. Need maybe mention if using complex 2D solid harmonics, formulas simpler; pseudo uses Cartesian. Let's think about performance of P2M direct to ancestors with atomics: each particle contributes to D nodes; D maybe 10; K 15; N=1e6 -> 150 million atomic doubles, too slow. Production: sort by leaf, compute leaf moments with block reduction, then bottom-up internal nodes with 4 children and center translation. Should mention strongly. In pseudo, note simple. Maybe we can write p2m bottom-up in pseudo to be more production-like. Let's see if we can implement bottom-up with center shift black-box. It may be cleaner: ``` def p2m(self): # leaf moments: one block per leaf, reduce particles for leaf in leaf_ids: M[leaf] = sum over particles in leaf of m * ... # internal bottom-up: for level reversed: for node in nodes_at_level(level): M[node] = 0 for child in children[node]: # translate child's moments about child center to node center M[node] += multipole_center_shift(M[child], center[child]-center[node], p) ``` Need `multipole_center_shift` formula. If child moments M_c,α about child center, parent center offset d = c_c - c_p. We want M_p,β = Σ_{α,γ: α+γ=β} M_c,α d^γ / γ! ? But sign? M_c,α = Σ m (-1)^{|α|} δ_child^α/α!. δ_parent = d + δ_child. M_p,β = Σ m (-1)^{|β|} (d+δ)^β/β!. Expand: Σ_α Σ_γ β!/(α!γ!) m (-1)^{|β|} d^γ δ^α. M_c,α = m (-1)^{|α|} δ^α/α!. So M_p,β = Σ_α,γ β!/(α!γ!) d^γ δ^α m (-1)^{|β|}. Express in M_c: M_c,α = m (-1)^{|α|} δ^α/α! => m δ^α = (-1)^{|α|} α! M_c,α. Then M_p,β = Σ_α,γ β!/(α!γ!) d^γ (-1)^{|α|+|β|} α! M_c,α = Σ_α,γ (-1)^{|β|} β!/γ! d^γ M_c,α. Since |α|+|γ|=|β|, (-1)^{|α|+|β|}=(-1)^{|2α? no|} = (-1)^{|β|} because 2|α| even. So M_p,β = Σ_{α+γ=β} (-1)^{|β|} β!/γ! M_c,α d^γ. Is that right? Test 1D: M_c0=m. d = child center offset. M_p1 = (-1)^1 1!/0! M_c0 d = -m d. But parent multipole M1 = Σ m (-1) (δ_parent)/1! = -m(d+δ) sum. If child moment M_c1 = -Σ m δ. M_p1 = -m d + M_c1? Formula gives α=1,γ=0: (-1)^1 1!/0! M_c1 = -M_c1? Wait β=1: α=0,γ=1 term: (-1)^1*1!/1!*M_c0*d = -m d. α=1,γ=0: (-1)^1*1!/0!*M_c1 = -M_c1. If M_c1=-Σmδ, -M_c1=Σmδ. Sum = Σm(δ-d). But expected M_p1 = -Σm(d+δ) = -md - Σmδ. Sign mismatch. Let's re-derive sign. Definition M_β = Σ m (-1)^{|β|} δ^β / β!. Child M_c,α = Σ m (-1)^{|α|} δ^α / α!. Parent δ_p = d + δ_c. M_p,β = Σ m (-1)^{|β|} (d+δ)^β / β! = Σ m (-1)^{|β|} Σ_{α+γ=β} β!/(α!γ!) d^γ δ^α / β! = Σ_α,γ (-1)^{|β|} / (α!γ!) d^γ δ^α m But m δ^α/α! = (-1)^{|α|} M_c,α. Thus M_p,β = Σ_α,γ (-1)^{|β|} (-1)^{|α|} / γ! d^γ M_c,α = Σ_α,γ (-1)^{|β|+|α|} d^γ/γ! M_c,α. Since |β|=|α|+|γ|, sign = (-1)^{2|α|+|γ|}=(-1)^{|γ|}. So M_p,β = Σ_{α+γ=β} (-1)^{|γ|} d^γ/γ! M_c,α. Test β=1: α=0,γ=1: (-1)^1 d M0 = -md. α=1,γ=0: M_c1 = -Σmδ. Sum = -md -Σmδ expected. Good. So center shift: `M_parent,β = Σ_{α+γ=β} (-1)^{|γ|} d^γ/γ! M_child,α`. Good. We can include this if bottom-up. But pseudo can just call `m2c(M_child, d)`. Maybe not necessary. Using bottom-up in pseudo may be more production and avoids atomics. Let's consider code: ``` def p2m(self): self.M.zero_() # leaf moments: parallel over leaves for leaf in self.tree.leaf_ids: s = torch.zeros(K) for i in range(tree.first[leaf], tree.last[leaf]): d = (x[i]-cx[leaf], y[i]-cy[leaf]) for k,alpha: s[k] += m[i]*... M[leaf]=s # bottom-up for level in range(max_depth-1, -1, -1): for node in tree.nodes_at_level(level): s = zeros(K) for child in children[node]: if child == -1: continue d = (cx[child]-cx[node], cy[child]-cy[node]) # shift child moments to node center for kb,beta in enumerate(ALPHAS): val=0 for ka,alpha in enumerate(ALPHAS): gamma = (beta[0]-alpha[0], beta[1]-alpha[1]) if gamma[0]>=0 and gamma[1]>=0: sign = -1 if (gamma[0]+gamma[1])%2 else 1 val += sign * M[child,ka] * monomial(d,gamma) / mf(gamma) s[kb] += val M[node]=s ``` This is more complex but still pseudo. It avoids atomic. Need `nodes_at_level` and `children`. Maybe too long. The direct ancestor pseudo is shorter. But for production explanation, mention bottom-up. Pseudo can use direct with note. I think direct is acceptable for pseudo code clarity. Need maybe include `far` list building in pseudo? We can leave as method stub. Good. Let's think about "Python pseudo code" maybe they expect high-level PyTorch/CuPy code, not CUDA C. Our code is pseudo with Python loops; okay. Could include a shorter alternative pseudo code using `cupy` vectorized for fixed grid? But FMM far list hard. Keep quadtree. Need maybe include "NVIDIA GPU" in pseudo code: `device='cuda'`. Good. Let's ensure final answer's equations are correct and not too many sign issues. We can say "one sign convention is; adjust if you define potential with negative sign". To avoid sign confusion, define acceleration directly and use formulas. Good. Let's write final answer in polished form. Potential final answer outline: ``` Assumption: Fast Multipole Method. 1. Approach - Use quadtree/FMM, 2D Newtonian kernel, softening. - Multipole/local expansion equations. - GPU implementation strategy. 2. Pseudo code ... code ... 3. Notes - true 2D log kernel, periodic, precision, validation. ``` Let's draft final text. Detailed explanation: **FMM math**: ``` K(r)=1/sqrt(r_x^2+r_y^2+eps^2) a_i = G Σ_j m_j ∇_{r_i} K(r_i - r_j) ``` Actually ∇K = (r_j-r_i)/r^3? K(r_i-r_j), gradient = -(r_i-r_j)/r^3 = (r_j-r_i)/r^3. Good. Multipole coefficients as above. M2P and M2L equations. Local polynomial evaluation. **Well separated**: `d_box` formula. **GPU**: - SoA, node arrays, M/L. - Kernels list with parallel mapping. - Tree rebuild. - Precision and softening. **Pseudo code** with comments. Let's ensure pseudo code's `inv_r_derivatives` returns D for all a,b up to max_order, but D dimensions `max_order+1` by `max_order+1`; accessing `D[alpha[0]+beta[0], alpha[1]+beta[1]]` where sum <= p. Good. Potential issue: In `m2l`, if `beta` total > p? ALPHAS up to p. For β=(p,0), α must 0. Good. Potential issue: In `l2p`, local coefficients `L` may have been accumulated for target leaf from far source nodes. If `L` includes contributions from source nodes that are also in near list due to far list error, double count. State far list excludes near. Potential issue: In `near_p2p`, if `near[b]` includes internal nodes not leaves? We use leaves. For adaptive, near list should include all leaf cells in neighborhood. Good. Potential issue: `tree.first` and `tree.last` for internal nodes? P2M direct uses leaf only. Near uses leaves. M2L uses far source nodes can be internal. M arrays for internal computed. Good. Maybe in pseudo code `self.tree.leaf_ids` property returns list of leaf node ids. `self.tree.range(node)` returns `(first,last)`. Good. Let's include `num_nodes` property. Maybe include `self.tree.nodes_at_level` not used if direct P2M. Not needed. Let's add comments in code: ``` # atomic_add is a placeholder; in CUDA use atomicAdd for float or double, # or better: sort particles by leaf and use block reductions. ``` Need mention `torch` is pseudo; real CUDA kernels use `float64_t`. Let's consider if we should include a small direct M2P version instead of M2L in pseudo to avoid far list standard bookkeeping. Maybe not. Could include a note after code: "The pseudo-code's `build_far_lists` is the main FMM bookkeeping step. In a production code, replace the simplified leaf-level M2L with the standard hierarchical M2L/L2L lists; the kernel structure remains the same." Good. Let's also include "If you use PyTorch, you can express leaf direct and L2P vectorized, but M2L/tree traversal is best in CUDA." Now, let's think about true 2D log kernel formulas. If user meant 2D gravity (infinite cylinders), force `a_i = 2πG Σ m_j (r_j-r_i)/(r^2+ε^2)` (since ∇ log r = r/r^2). FMM can use `K(r)=log(r)`; multipole expansion in 2D has logarithmic local/multipole terms. The pseudo code's `inv_r_derivatives` would be `log_derivatives`. Mention. Potential answer length: likely long. Desired oververbosity 9, okay. Let's write final now. Ensure no analysis in final. Use Markdown. Include code block. Double-check equations: - `M_A,α = Σ m (-1)^|α| δ^α/α!`. - M2P: `a += G Σ M_α ∇ ∂^α K(R)`. In code, `D[α+e]` is ∂^{α+e} K. Good. - M2L: `C_β += Σ_{α+β≤p} M_α ∂^{α+β} K(D) / β!`. In code, `Dtab[α+β]/mf(beta)`. Good. - L2P: `a += G Σ C_β ∇ η^β`. In code, `grad_monomial`. Good. Test single source: M0=m, M1=-mδ, ... C1 = M0 ∂K + M1 ∂²K/1!; L2P C1 constant acceleration. Good. One subtlety: In M2L, `D = c_B - c_A`. `K(D) = 1/|c_B-c_A|`. But source expansion originally uses R = r_i - c_A = D + η. Taylor `K(D+η)` expansion uses derivatives at D with positive η. Our local formula uses `∂^{α+β} K(D) η^β/β!` no (-1)^β. Good. Test 1D D>0, K(D+η)=1/D - η/D² + η²/D³... C1 from M0 ∂K = -m/D²; L2P C1 constant -m/D². C2 from M0 ∂²K/2 = (2/D³)/2=1/D³; L2P 2η/D³; expansion acceleration ∂η K = -1/D² + 2η/D³... good. No sign issue. Multipole M sign: Source expansion `K(R) = Σ M_α ∂^α K(R)` with M_α = Σ m (-1)^|α| δ^α/α!`. Test source at δ, target R (relative A): K=1/|R-δ| = 1/R + δ/R² + δ²/R³... M0=m, M1=-mδ. M2P potential? K approx M0 ∂0 K + M1 ∂1 K = m/R + (-mδ)*(-1/R²)=m/R + mδ/R². Good. Acceleration ∇: G [ -m/R² + 2mδ/R³? Wait derivative of mδ/R² w.r.t R = -2mδ/R³. But actual acceleration ∂R 1/(R-δ) = -1/(R-δ)² = -1/R² + 2δ/R³... Our M1 term: M1 ∂ (∂K) = M1 ∂²K = (-mδ)*(2/R³) = -2mδ/R³. That gives -1/R² -2δ/R³. But expected for target R relative A, source at δ: distance R-δ, acceleration ∂R 1/(R-δ) = -1/(R-δ)² ≈ -1/R² + 2δ/R³? Let's compute: (R-δ)^-2 = R^-2 (1-δ/R)^-2 = R^-2 (1+2δ/R+3δ²/R²) = 1/R²+2δ/R³. With negative sign: -1/R² -2δ/R³. Yes expected negative. Good. M1 term negative. Good. M2L with D = c_B-c_A, source δ. C1 = M0 ∂K(D) + M1 ∂²K(D) = -m/D² -2mδ/D³. L2P constant. Actual acceleration at target near B: distance D+η-δ, acceleration ∂η 1/(D+η-δ) = -1/(D+η-δ)² ≈ -1/D² -2δ/D³ + ... Good. All consistent. Need maybe include `G` in M2L? We store without G, apply in L2P/near. Good. Let's consider if local coefficients should be accumulated over source nodes `a` in far list. If far list includes source nodes at different depths, M_A includes all particles in node. Good. Now, let's write final with enough detail. Potential final code: Use `self.eps2 = self.eps*self.eps`? In direct, use `eps*eps`. In `inv_r_derivatives`, pass eps. Good. Add `self.theta`, `self.p` stored. Add `self.step_count = 0`. Add `domain` maybe from initial positions with margin. Let's write code elegantly. One more thought: In P2M direct to ancestors, if tree has internal nodes with children, a particle contributes to all ancestors. But if a node has many particles, atomic adds. In pseudo okay. But M for internal node computed from all particles in subtree; good. M for leaf also. M2L can use internal node M. Near direct uses leaf particles. Good. Far list for target leaf may include internal nodes that are ancestors? Exclude. If an internal node is not well-separated but its child might be? In standard traversal, if not well-separated, descend. Our far list should be built by traversal: for target leaf, traverse source tree; if node well-separated, add to far; else if leaf, add to near; else recurse. This ensures all source particles are either in near leaves or far multipoles. Good. We can state this in explanation and pseudo `build_far_lists` comment. That resolves correctness. The near list can be built as leaves reached by recursion when not well-separated. For a target leaf, near leaves are those in the subtree of not-far nodes. It may not be just 3x3 at same depth; adaptive near includes finer cells. Good. So in explanation: "For each target leaf, do a tree walk: if source node is well-separated, record it for M2L; otherwise if it is a leaf, record it for P2P; otherwise recurse." This is correct and simple. Complexity? For each target leaf, traversal could be O(N) in worst if many near? But near bounded by leaf capacity and θ; far stops early. Standard O(N log N). Good. This far list per target leaf may be large but okay. It avoids L2L double count. It is essentially Barnes-Hut-like M2L to leaf local expansion. Is it FMM? It uses local expansion at leaf only, not hierarchical L2L. But it's a multipole-accelerated method. To call it FMM, we can say "A full production FMM adds local expansion propagation (L2L) and multi-level M2L to reduce far-list size; the pseudo-code uses the simpler leaf-local version for readability." Good. Maybe in explanation, present standard FMM stages, then in pseudo code use leaf-local M2L. That's honest. Let's phrase: "The mathematically complete FMM is P2M, M2L, L2L, L2P, near P2P. The pseudo-code below shows the essential kernels; for brevity it builds a leaf-level local expansion from a precomputed far list (equivalent to doing M2L directly for leaves). A production NVIDIA implementation replaces that direct leaf far list with the standard multi-level M2L/L2L list to reduce constant factors." Good. Need maybe include L2L in explanation but not code? Could include optional code snippet for L2L? Maybe too long. We can mention. Let's include a short L2L formula in explanation: If local acceleration polynomial `a(η)=Σ C_β ∇η^β` is shifted by `d`, then `C'_γ = Σ_{β≥γ} C_β (β!/(γ!(β-γ)!)) d^{β-γ}`? Wait identity earlier: ∇(d+η)^β = Σ_{γ≤β} β!/(γ!(β-γ)!) d^{β-γ} ∇η^γ. Coefficient for C_β to C'_γ is multinomial(β,γ) d^{β-γ}. No extra? For β=2, γ=1: 2 d, correct. For β=3, γ=1: 3 d^2? ∇(d+η)^3 = 3(d+η)^2 = 3d^2 +6dη+3η². ∇η^1=1, ∇η^2=2η, ∇η^3=3η². Coeff: γ=1:3d²; γ=2:3d? because 3d * ∇η² =6dη; γ=3:1*∇η³=3η². Multinomial(3,1)=3, (3,2)=3, (3,3)=1. Good. So L2L formula: `C_child,γ += C_parent,β * multinomial(β,γ) * d^{β-γ}`. We can include in explanation as standard. But in our local acceleration coefficients, C_β is coefficient for ∇η^β. The shift identity holds. Good. Now, let's write final. Need ensure final answer says "Python pseudo code" not necessarily executable. Good. Let's produce final. Below I assume that “fast multiple method” means **Fast Multipole Method (FMM)**. I will design a 2D N-body gravity simulator for NVIDIA GPUs using a quadtree FMM. The pseudo-code is Python-like, but the loops indicate CUDA thread/block work, not literal Python loops. --- # 1. Physical model We have \(N\) point particles in the \(x,y\) plane: \[ \mathbf r_i=(x_i,y_i), \qquad m_i \] For standard Newtonian gravity restricted to 2D positions, use the 3D inverse-distance kernel in the plane: \[ K(\mathbf r)=\frac{1}{\sqrt{r_x^2+r_y^2+\epsilon^2}} \] The acceleration of particle \(i\) is \[ \mathbf a_i = G \sum_{j\ne i} m_j \nabla_{\mathbf r_i} K(\mathbf r_i-\mathbf r_j) = G \sum_{j\ne i} m_j \frac{\mathbf r_j-\mathbf r_i} {\left(|\mathbf r_i-\mathbf r_j|^2+\epsilon^2\right)^{3/2}} \] where \(G\) is the gravitational constant and \(\epsilon\) is a softening length. > If by “2D gravity” you mean a strictly 2D Poisson/logarithmic kernel, replace \(K(r)=1/r\) with \(K(r)=\ln r\). The FMM structure is the same, but the derivative tables and translation kernels change. --- # 2. Fast Multipole Method idea The direct method is \(O(N^2)\). FMM separates interactions into: 1. **Near field**: particle-to-particle, direct. 2. **Far field**: cluster-to-particle, multipole expansion. A spatial quadtree partitions the domain. Each node/cluster has: - center \(\mathbf c_A\) - half-size \(h_A\) - multipole coefficients \(M_{A,\alpha}\) - local expansion coefficients \(L_{B,\beta}\) Use multi-indices: \[ \alpha=(\alpha_x,\alpha_y), \qquad |\alpha|=\alpha_x+\alpha_y \] For a source cluster \(A\), define multipole coefficients up to degree \(p\): \[ M_{A,\alpha} = \sum_{j\in A} m_j \frac{(-1)^{|\alpha|}}{\alpha!} (\mathbf r_j-\mathbf c_A)^\alpha \] where \[ \alpha! = \alpha_x!\alpha_y! \] and \[ (\mathbf r_j-\mathbf c_A)^\alpha = (r_{jx}-c_{Ax})^{\alpha_x} (r_{jy}-c_{Ay})^{\alpha_y} \] --- ## 2.1 Multipole-to-particle, M2P For a target particle \(i\), let \[ \mathbf R = \mathbf r_i-\mathbf c_A \] The far-field acceleration from cluster \(A\) is approximated by \[ \mathbf a_i^{(A)} \approx G \sum_{|\alpha|\le p} M_{A,\alpha} \nabla \partial^\alpha K(\mathbf R) \] where \[ \partial^\alpha K = \partial_x^{\alpha_x} \partial_y^{\alpha_y} K \] In practice, this means computing derivatives of \[ K(\mathbf R)=\frac{1}{\sqrt{R_x^2+R_y^2+\epsilon^2}} \] up to order \(p+1\). --- ## 2.2 Multipole-to-local, M2L For a target cluster \(B\) with center \(\mathbf c_B\), define \[ \mathbf D = \mathbf c_B-\mathbf c_A \] We build local acceleration coefficients \(L_{B,\beta}\) such that particles inside \(B\) can evaluate the far field cheaply. One convention is: \[ L_{B,\beta} \mathrel{+}= \sum_{|\alpha|+|\beta|\le p} \frac{ M_{A,\alpha} }{ \beta! } \partial^{\alpha+\beta}K(\mathbf D) \] Then for a particle \(i\in B\), let \[ \boldsymbol\eta = \mathbf r_i-\mathbf c_B \] The local expansion contributes \[ \mathbf a_i^{(B)} = G \sum_{|\beta|\le p,\ \beta\ne 0} L_{B,\beta} \nabla \boldsymbol\eta^\beta \] This is a small polynomial evaluation, very GPU-friendly. --- ## 2.3 Local-to-local, L2L In a full hierarchical FMM, local expansions are propagated from parent nodes to child nodes. If a parent has local coefficients \(L_{\text{parent},\beta}\), and the child center is displaced by \(\mathbf d\), then \[ L_{\text{child},\gamma} \mathrel{+}= \sum_{\beta\ge \gamma} L_{\text{parent},\beta} \frac{\beta!}{\gamma!(\beta-\gamma)!} \mathbf d^{\beta-\gamma} \] This is just a polynomial shift. --- # 3. Well-separated criterion A source node \(A\) is used as a multipole approximation for target node/particle \(B\) if they are **well separated**. A practical criterion: Let the minimum box distance between \(A\) and \(B\) be \(d_{\text{box}}\). Use multipole if \[ d_{\text{box}} > \theta (h_A+h_B) \] where: - \(\theta\) is the opening-angle parameter, often \(0.5\) to \(1.0\) - smaller \(\theta\) gives more accuracy but more work - larger multipole degree \(p\) gives more accuracy but larger translations For a quadtree, a common practical rule is: - If the source node is in the 3×3 neighborhood of the target cell, use direct P2P. - Otherwise, if the well-separated criterion holds, use multipole. - Otherwise, recurse into children. --- # 4. GPU implementation strategy ## 4.1 Data layout Use structure-of-arrays for particles: ```python x, y, vx, vy, ax, ay, m ``` On GPU as `float32` or `float64`. For high accuracy, use `float64` for moments and local expansions; direct near-field can often be `float64` as well. Tree node arrays: ```python cx, cy, h, parent, child[4], first_particle, last_particle ``` FMM arrays: ```python M[node, k] # multipole moments L[node, k] # local expansion coefficients ``` where \[ k = 0,\dots,K-1 \] and \[ K = \frac{(p+1)(p+2)}{2} \] For example: - \(p=4 \Rightarrow K=15\) - \(p=6 \Rightarrow K=28\) --- ## 4.2 Main kernels ### `P2M`: particles to multipole Compute multipole moments for every tree node. Conceptually: - each particle contributes to its leaf and ancestors - in production, sort particles by leaf, use block reductions, then bottom-up accumulation --- ### `M2L`: multipole to local For each target node/leaf, loop over well-separated source nodes and translate their multipole coefficients into local coefficients. This kernel is highly parallel over target nodes. Inside each block: - load source multipole coefficients into shared memory - evaluate translation kernels - accumulate local coefficients For small \(p\), the translation matrices are small. For production performance: - precompute translation tables when possible - unroll loops for \(p=4,6,8\) - use double precision or mixed precision with error correction --- ### `L2L`: local to local Top-down propagation of local expansions from parent to children. Parallel over nodes. --- ### `L2P`: local to particle One thread per particle. Load the local expansion of the particle’s leaf into registers/shared memory and evaluate the polynomial. --- ### `P2P`: near-field particle-particle One CUDA block per leaf cell. Load the particles in the leaf and neighboring leaves into shared memory. For each particle in the leaf: - loop over nearby particles - compute direct softened acceleration If leaf capacity is `leaf_max`, the near-field cost is roughly \[ O(N \cdot \text{leaf\_max}) \] --- ## 4.3 Tree rebuilding Because particles move, the quadtree must be updated. Options: 1. Rebuild every timestep. 2. Rebuild every \(K\) timesteps. 3. Rebuild only when a cell becomes overloaded or particles cross boundaries. For GPU implementation: - compute max-depth cell ID for each particle - radix sort particles by cell ID - count particles per cell - build occupied quadtree nodes - store leaf assignment and particle ranges Rebuilding can be done on GPU, but many simulations rebuild on CPU and copy tree metadata to GPU every few steps. --- # 5. Python pseudo-code This is conceptual pseudo-code for an NVIDIA GPU implementation. The loops are meant to show the parallel work distribution. A real implementation would use CUDA C++, CuPy, PyTorch custom ops, CUTLASS/CuTe, or similar. ```python # ===================================================================== # Python pseudo-code for a 2D N-body FMM on NVIDIA GPU # # Not production code. # Loops indicate CUDA thread/block work. # Use float64 for accuracy, float32 if acceptable. # ===================================================================== import torch from math import factorial, hypot G = 1.0 # gravitational constant; set to 1 for pseudo-code # --------------------------------------------------------------------- # Multi-index helpers # --------------------------------------------------------------------- def multi_indices(p): """ Return all 2D multi-indices alpha=(a,b) with a+b <= p. """ return [ (a, b) for a in range(p + 1) for b in range(p + 1 - a) ] def multi_factorial(alpha): return factorial(alpha[0]) * factorial(alpha[1]) def monomial(v, alpha): return v[0] ** alpha[0] * v[1] ** alpha[1] def grad_monomial(v, alpha): """ gradient of v^alpha """ ax = 0.0 ay = 0.0 if alpha[0] > 0: ax = ( alpha[0] * v[0] ** (alpha[0] - 1) * v[1] ** alpha[1] ) if alpha[1] > 0: ay = ( alpha[1] * v[0] ** alpha[0] * v[1] ** (alpha[1] - 1) ) return ax, ay # --------------------------------------------------------------------- # Kernel derivative table # --------------------------------------------------------------------- def inv_r_derivatives(R, max_order, eps): """ Return D[a,b] = partial_x^a partial_y^b K(R) where K(R) = 1 / sqrt(Rx^2 + Ry^2 + eps^2) for 0 <= a+b <= max_order. In production: - precompute symbolic polynomial formulas for p=4,6,8 - or generate them offline with sympy - or use a recursive polynomial table Useful recurrence: Let r2 = Rx^2 + Ry^2 + eps^2 Let D[a,b] = r ** (-(2*(a+b)+1)) * Q[a,b](Rx,Ry) Q[0,0] = 1 Q[a,b] = r2 * dQdx(a-1,b) - (2*(a+b)-1) * Rx * Q[a-1,b] Q[a,b] = r2 * dQdy(a,b-1) - (2*(a+b)-1) * Ry * Q[a,b-1] For small p, unroll these formulas in CUDA. """ D = [ [0.0] * (max_order + 1) for _ in range(max_order + 1) ] # Pseudo: fill D using precomputed formulas. # A real CUDA implementation would store this in registers/shared memory. return D # --------------------------------------------------------------------- # Quadtree # --------------------------------------------------------------------- class QuadTree: def __init__(self, domain, max_depth, leaf_max): """ domain = (xmin, ymin, xmax, ymax) """ self.max_depth = max_depth self.leaf_max = leaf_max # Pseudo: node arrays on GPU. # In production, allocate to max_nodes and track num_nodes. self.cx = torch.zeros(1, dtype=torch.float64, device="cuda") self.cy = torch.zeros(1, dtype=torch.float64, device="cuda") self.h = torch.zeros(1, dtype=torch.float64, device="cuda") self.parent = torch.zeros(1, dtype=torch.int32, device="cuda") self.child = torch.zeros((1, 4), dtype=torch.int32, device="cuda") # particle ranges for each node self.first = torch.zeros(1, dtype=torch.int32, device="cuda") self.last = torch.zeros(1, dtype=torch.int32, device="cuda") # leaf assignment for each particle self.leaf = torch.zeros(1, dtype=torch.int32, device="cuda") self.num_nodes = 1 def assign(self, x, y): """ Build or rebuild quadtree from current particle positions. Production GPU approach: 1. compute max-depth grid cell for each particle 2. radix sort particles by cell 3. count particles per cell 4. build occupied quadtree nodes 5. subdivide cells with count > leaf_max 6. store parent/child/first/last/leaf arrays For pseudo-code, assume this works and returns node arrays. """ pass def leaf_ids(self): """ Return list/array of leaf node IDs. """ pass def range(self, node): """ Return (first_particle, last_particle) for node. """ return self.first[node], self.last[node] def well_separated(self, a, b, theta): """ Conservative box-distance well-separated test. """ dx = max( 0.0, abs(self.cx[a] - self.cx[b]) - self.h[a] - self.h[b] ) dy = max( 0.0, abs(self.cy[a] - self.cy[b]) - self.h[a] - self.h[b] ) d_box = hypot(dx, dy) return d_box > theta * (self.h[a] + self.h[b]) def build_near_lists(self): """ For each target leaf, build near leaf list. Standard method: For target leaf B: traverse tree from root: if source node A is well-separated from B: do not include its descendants in near list else if A is a leaf: add A to near[B] else: recurse into children of A In a fixed-depth quadtree, this is approximately the 3x3 neighborhood. In an adaptive quadtree, it includes finer neighboring leaves. """ pass def build_far_lists(self, theta): """ For each target leaf, build far source list for M2L. Standard method: For target leaf B: traverse tree from root: if source node A is well-separated from B: add A to far[B] else: recurse into children of A The resulting far list contains internal nodes and/or leaves whose particles are far enough to be represented by multipole. """ pass # --------------------------------------------------------------------- # FMM N-body simulator # --------------------------------------------------------------------- class NBody2D_FMM: def __init__( self, x, y, m, domain, max_depth=8, p=4, theta=0.7, eps=1e-3, leaf_max=64 ): """ x, y, m: initial particle positions and masses domain: (xmin, ymin, xmax, ymax) p: multipole/local expansion degree theta: FMM opening angle parameter eps: softening leaf_max: maximum particles per leaf cell """ # Particle arrays on GPU. self.x = torch.tensor(x, dtype=torch.float64, device="cuda") self.y = torch.tensor(y, dtype=torch.float64, device="cuda") self.m = torch.tensor(m, dtype=torch.float64, device="cuda") self.N = self.x.numel() self.vx = torch.zeros_like(self.x) self.vy = torch.zeros_like(self.y) self.ax = torch.zeros_like(self.x) self.ay = torch.zeros_like(self.y) self.p = p self.theta = theta self.eps = eps self.eps2 = eps * eps self.leaf_max = leaf_max # Tree. self.tree = QuadTree(domain, max_depth, leaf_max) self.tree.assign(self.x, self.y) # FMM index table. self.ALPHAS = multi_indices(p) self.K = len(self.ALPHAS) # Node FMM arrays. self.M = torch.zeros( (self.tree.num_nodes, self.K), dtype=torch.float64, device="cuda" ) self.L = torch.zeros_like(self.M) # FMM neighbor lists. self.far = self.tree.build_far_lists(self.theta) self.near = self.tree.build_near_lists() self.step_count = 0 # ----------------------------------------------------------------- # P2M: particles to multipole # ----------------------------------------------------------------- def p2m(self): """ Compute multipole moments for all tree nodes. Pseudo simple version: Each particle contributes to its leaf and all ancestors. Production version: 1. sort particles by leaf 2. compute leaf moments with block reductions 3. bottom-up compute internal moments with center shifts The pseudo version is easier to read but uses atomicAdd in a real CUDA kernel. """ self.M.zero_() # CUDA kernel pseudo: # one thread per particle # depth is small, so ancestor loop is short # use atomicAdd(double) or better: block reductions for i in range(self.N): node = self.tree.leaf[i] while node >= 0: dx = self.x[i] - self.tree.cx[node] dy = self.y[i] - self.tree.cy[node] for k, alpha in enumerate(self.ALPHAS): a, b = alpha sign = -1.0 if (a + b) % 2 == 1 else 1.0 val = ( sign * self.m[i] * (dx ** a) * (dy ** b) / multi_factorial(alpha) ) # CUDA pseudo: atomicAdd(self.M[node, k], val) self.M[node, k] = self.M[node, k] + val node = self.tree.parent[node] # ----------------------------------------------------------------- # M2L: multipole to local # ----------------------------------------------------------------- def m2l(self): """ Compute local expansion coefficients for target leaves from well-separated source nodes. CUDA kernel pseudo: one block per target leaf block loops over far[source] list source M coefficients loaded into shared memory translation accumulated into L[target] """ self.L.zero_() for b in self.tree.leaf_ids(): # b is target leaf for a in self.far[b]: # a is well-separated source node D = ( self.tree.cx[b] - self.tree.cx[a], self.tree.cy[b] - self.tree.cy[a] ) # Derivatives of K(D) up to order p. Dtab = inv_r_derivatives(D, self.p, self.eps) for kb, beta in enumerate(self.ALPHAS): # beta = local polynomial index if beta == (0, 0): # constant local acceleration term is not used continue acc = 0.0 for ka, alpha in enumerate(self.ALPHAS): if ( alpha[0] + alpha[1] + beta[0] + beta[1] ) <= self.p: acc += ( self.M[a, ka] * Dtab[ alpha[0] + beta[0], alpha[1] + beta[1] ] / multi_factorial(beta) ) # CUDA pseudo: atomicAdd(self.L[b, kb], acc) self.L[b, kb] = self.L[b, kb] + acc # ----------------------------------------------------------------- # Optional L2L: local to local # ----------------------------------------------------------------- def l2l(self): """ Top-down local-to-local translation. This is included conceptually for completeness. The simplified pseudo-code above builds leaf local expansions directly from far lists; a production FMM often uses L2L to reduce the number of M2L operations. """ # Top-down over tree levels. for level in range(self.tree.max_depth): # CUDA pseudo: one thread/block per node at this level for node in self.tree.nodes_at_level(level): for c in range(4): child = self.tree.child[node, c] if child < 0: continue d = ( self.tree.cx[child] - self.tree.cx[node], self.tree.cy[child] - self.tree.cy[node] ) # Shift local polynomial from parent to child. for kb, beta in enumerate(self.ALPHAS): for kg, gamma in enumerate(self.ALPHAS): # gamma <= beta componentwise if ( gamma[0] <= beta[0] and gamma[1] <= beta[1] ): # multinomial(beta, gamma) coef = ( factorial(beta[0]) / ( factorial(gamma[0]) * factorial(beta[0] - gamma[0]) ) ) * ( factorial(beta[1]) / ( factorial(gamma[1]) * factorial(beta[1] - gamma[1]) ) ) dcoef = monomial( d, ( beta[0] - gamma[0], beta[1] - gamma[1] ) ) self.L[child, kg] += ( self.L[node, kb] * coef * dcoef ) # ----------------------------------------------------------------- # L2P: local to particle # ----------------------------------------------------------------- def l2p(self): """ Evaluate local expansions at particle positions. CUDA kernel pseudo: one thread per particle load L[leaf] into registers/shared memory evaluate small polynomial """ for i in range(self.N): b = self.tree.leaf[i] eta = ( self.x[i] - self.tree.cx[b], self.y[i] - self.tree.cy[b] ) ax = 0.0 ay = 0.0 for k, beta in enumerate(self.ALPHAS): if beta == (0, 0): continue gx, gy = grad_monomial(eta, beta) ax += self.L[b, k] * gx ay += self.L[b, k] * gy # Add far-field local expansion contribution. self.ax[i] += G * ax self.ay[i] += G * ay # ----------------------------------------------------------------- # P2P: near-field direct interaction # ----------------------------------------------------------------- def near_p2p(self): """ Direct particle-particle near field. CUDA kernel pseudo: one block per leaf cell load this leaf and neighboring leaves into shared memory one thread per particle in the leaf loop over nearby particles """ for b in self.tree.leaf_ids(): nb_list = self.near[b] # Production: # load x, y, m for all particles in nb_list # into shared memory before the particle loop. i0, i1 = self.tree.range(b) for i in range(i0, i1): ax = 0.0 ay = 0.0 for nb in nb_list: j0, j1 = self.tree.range(nb) for j in range(j0, j1): if i == j: continue dx = self.x[j] - self.x[i] dy = self.y[j] - self.y[i] r2 = dx * dx + dy * dy + self.eps2 inv_r3 = r2 ** (-1.5) ax += self.m[j] * dx * inv_r3 ay += self.m[j] * dy * inv_r3 self.ax[i] += G * ax self.ay[i] += G * ay # ----------------------------------------------------------------- # Full FMM acceleration # ----------------------------------------------------------------- def compute_accelerations(self): """ Compute acceleration for all particles using FMM. """ self.ax.zero_() self.ay.zero_() self.p2m() self.m2l() # In a full hierarchical FMM: # self.l2l() # before L2P. self.l2p() self.near_p2p() # ----------------------------------------------------------------- # Velocity Verlet timestep # ----------------------------------------------------------------- def step(self, dt, rebuild_every=10): """ One velocity Verlet step. For production: - use CUDA streams - avoid host synchronization - copy new tree only when rebuilt - use double precision for FMM moments/local expansions """ a0x = self.ax.clone() a0y = self.ay.clone() # Position update. self.x += self.vx * dt + 0.5 * a0x * dt * dt self.y += self.vy * dt + 0.5 * a0y * dt * dt # Rebuild tree occasionally. if self.step_count % rebuild_every == 0: self.tree.assign(self.x, self.y) self.far = self.tree.build_far_lists(self.theta) self.near = self.tree.build_near_lists() # If tree size changed, reallocate M and L. # Pseudo-code ignores reallocation. # New acceleration. self.compute_accelerations() # Velocity update. self.vx += 0.5 * (a0x + self.ax) * dt self.vy += 0.5 * (a0y + self.ay) * dt self.step_count += 1 # ===================================================================== # Example usage # ===================================================================== if __name__ == "__main__": # N = 100000 # x = ... # y = ... # m = ... # domain = (xmin, ymin, xmax, ymax) # sim = NBody2D_FMM( # x=x, # y=y, # m=m, # domain=domain, # max_depth=8, # p=4, # theta=0.7, # eps=1e-3, # leaf_max=64 # ) # sim.compute_accelerations() # for step in range(10000): # sim.step(dt=1e-4) pass ``` --- # 6. Practical choices ## Multipole degree \(p\) Good starting values: | Accuracy | Suggested \(p\) | |---|---:| | Fast / low accuracy | 3 or 4 | | General purpose | 5 or 6 | | High accuracy | 7 or 8 | Larger \(p\) increases: - FMM array size - translation cost - local evaluation cost but reduces truncation error. --- ## Opening angle \(\theta\) Smaller \(\theta\) is more accurate because far clusters are treated more conservatively. Typical values: \[ \theta = 0.5,\ 0.7,\ 1.0 \] A good test is: 1. Run with \(p=4, \theta=0.7\) 2. Run with \(p=6, \theta=0.5\) 3. Compare accelerations or energy error If they agree to required tolerance, the first may be acceptable. --- ## Leaf size `leaf_max` controls near-field cost. Typical values: ```text 32, 64, 128 ``` Small leaf size: - cheaper near field - more tree nodes - larger FMM overhead Large leaf size: - fewer nodes - more expensive near field - potentially larger local error if cell is large --- # 7. Validation For a new implementation, always validate against a direct O(N²) simulation for small \(N\). Check: 1. **Force/acceleration error** \[ \left\|\mathbf a_{\text{FMM}} - \mathbf a_{\text{direct}}\right\| \] 2. **Energy drift** \[ E = \sum_i \frac{1}{2}m_i|\mathbf v_i|^2 + U \] For an isolated system, total energy should drift slowly. 3. **Convergence with \(p\)** Error should decrease as \(p\) increases. 4. **Convergence with \(\theta\)** Smaller \(\theta\) should reduce error. 5. **Convergence with leaf size** Smaller leaf size should reduce near-field truncation error. --- # 8. Important implementation notes ## Double counting The FMM list construction must ensure each source contribution is used exactly once: - near sources: direct P2P - far sources: M2L/L2P - not both The tree-walk rule ```text if well_separated(source, target): use multipole else: descend ``` ensures this when applied consistently. --- ## Self-interaction Always skip \(i=j\) in direct P2P. The FMM far field should never include the particle’s own cell or ancestor cells in a way that includes the particle itself. --- ## Softening Use the softened kernel everywhere: \[ K(r)=\frac{1}{\sqrt{r^2+\epsilon^2}} \] This prevents singularities and makes derivatives stable. --- ## Periodic boundaries The design above assumes an open boundary or a fixed non-periodic domain. For periodic boundaries: - either use a periodic FMM - or use Ewald summation - or use a non-periodic FMM with image particles, which is expensive Periodic N-body FMM is a separate, more complex problem. --- ## NVIDIA GPU performance tips 1. Use structure-of-arrays. 2. Store FMM coefficients as 2D arrays `[node, K]`. 3. Use `float64` for moments/local expansions if accuracy matters. 4. Precompute translation tables when possible. 5. Unroll FMM kernels for fixed \(p\). 6. Use shared memory for near-field neighbor particles. 7. Use radix sort for particle-to-leaf assignment. 8. Avoid host-device synchronization inside the timestep. 9. Use CUDA graphs if the sequence of kernels is fixed and tree metadata is static for several steps. 10. Profile with Nsight Compute to check occupancy, memory bandwidth, and atomic contention.

FMM N-Body Gravity Sim