FMM M2M-Übersetzungsdiskrepanz (Greengards Kurzkurs) mit SciPy Spherical Harmonics ⇐ Python

[Anonymous]

Ich wende die Fast Multipole Method (FMM) für Coulomb-Wechselwirkungen an und folge „A Short Course on the Fast Multipole Method“ von Leslie Greengard: http://math.nyu.edu/faculty/greengar/sh ... se_fmm.pdf.
Mein Ziel ist es, das Potenzial für ein System von Punktladungen genau zu berechnen.
Arbeitskomponenten (P2M und M2P):
Meine Implementierung des Particle-to-Multipole (P2M)-Schritts und die anschließende Bewertung des Potenzials aus einer Multipolentwicklung (Multipol-to-Particle, M2P) scheint korrekt zu sein.

• Das Potenzial basiert auf Gleichung (5-15) aus dem Kurzkurs: Phi(P) = sum_n=0^inf sum_m=-n^n (M_n^m / r^(n+1)) * Y_n^m(theta, phi).
• Ich verwende scipy.special.sph_harm(m, n, phi, theta) für sphärische Harmonische, was orthonormales Y_n^m(theta, phi) zurückgibt.
• Meine P2M-Funktion berechnet Multipolmomente M_n^m (nennen wir sie M_P2M) als Summe q_i * r_i^n * conj(Y_n^m(theta_i, phi_i)).
• Meine Auswertungsfunktion verwendet diese M_P2M-Momente und wendet den Faktor (4 * pi / (2n + 1)) korrekt an, im Einklang mit orthonormalen Harmonischen.

Dies wird durch die folgende Ausgabe bestätigt:

Code: Select all

Potential from child multipole: 3.275435467795e+00
Potential from direct sum: 3.275435467795e+00

Wie Sie sehen können, stimmt das aus der Multipolexpansion (P2M) einer untergeordneten Zelle berechnete Potenzial sehr gut mit der direkten Summenberechnung überein.
Das Problem (M2M-Übersetzung):
Ich stoße jedoch auf Probleme mit der Multipol-zu-Multipol-Übersetzung (M2M), die durch Theorem 5.3 (Gleichung 5-20) in definiert ist Papier. Nach der Übertragung der Multipolmomente des Kindes auf das Zentrum der Elternzelle ist das resultierende Potenzial deutlich anders:

Code: Select all

Potential from parent multipole: 1.123673185858e+00
Potential from direct sum: 3.275435467795e+00

Dies weist auf einen Fehler in der M2M-Übersetzung hin. Zur Verifizierung berechne ich auch die Multipolmomente der Elternzelle direkt aus den Quellpartikeln (d. h. ich behandle die Elternzelle als neues Kind und führe P2M für ihre Quellen durch). Die Werte von parent_multipole_M2M (berechnet über M2M-Übersetzung) sind nicht gleich den Werten von parent_multipole_direct (berechnet direkt über P2M auf denselben Quellen, aber zentriert auf das übergeordnete Element). Ich glaube, der Fehler liegt in meiner Funktion M2M_theorem53.
Ich habe den gesamten Code überprüft, kann den Fehler jedoch nicht finden. Kann mir jemand helfen? Hier ist mein Code:

Code: Select all

from scipy.special import sph_harm, factorial
import numpy as np

def A_nm(n, m):
"""Compute A_n^m coefficient as defined in Greengard's paper."""
return (-1)**n / np.sqrt(factorial(n - abs(m)) * factorial(n + abs(m)))

def M2M_theorem53(child_multipole, child_center, parent_center, p):
"""
Perform M2M translation based on Theorem 5.3 (Greengard, 1997).
Translates multipole expansion from child center to parent center (origin).

Parameters:
child_multipole: np.array of shape (p+1, 2p+1), with child moments O_n^m at [n, m+n]
child_center: Point object
parent_center: Point object
p: order of expansion

Returns:
parent_multipole: np.array of shape (p+1, 2p+1), with translated M_j^k
"""
parent_multipole = np.zeros((p + 1, 2 * p + 1), dtype=complex)

# Vector from parent to child center
dx = child_center.x - parent_center.x
dy = child_center.y - parent_center.y
dz = child_center.z - parent_center.z

rho = np.sqrt(dx**2 + dy**2 + dz**2)
if rho < 1e-14:
return np.copy(child_multipole)  # No translation needed

alpha = np.arctan2(dy, dx)
beta = np.arccos(dz / rho)

for j in range(p + 1):
for k in range(-j, j + 1):
acc = 0.0j
A_jk = A_nm(j, k)
for n in range(j + 1):
for m in range(-n, n + 1):
j_minus_n = j - n
k_minus_m = k - m

if abs(k_minus_m) > j_minus_n:
continue  # Spherical harmonic is 0 outside its domain

A_nm_val = A_nm(n, m)
A_jn_km = A_nm(j_minus_n, k_minus_m)

# Complex exponential factor: i^{|k| - |m| - |k - m|}
phase = 1j ** (abs(k) - abs(m) - abs(k - m))

# Y_n^{-m}(alpha, beta)
Y = sph_harm(-m, n, alpha, beta)

term = (child_multipole[j_minus_n, k_minus_m + j_minus_n] *
phase * A_nm_val * A_jn_km * (rho ** n) * Y / A_jk)
acc += term
parent_multipole[j, k + j] = acc

return parent_multipole

import numpy as np
from scipy.special import sph_harm
# --- Re-using your Point and Particle classes ---
class Point():
"""The class for a point."""
def __init__(self, coords=[], domain=1.0):
if coords:
assert len(coords) == 3, "the size of coords should be 3."
self.x = coords[0]
self.y = coords[1]
self.z = coords[2]
else:
self.x = domain * np.random.random()
self.y = domain * np.random.random()
self.z = domain * np.random.random()

def distance(self, other):
return np.sqrt((self.x-other.x)**2 + (self.y-other.y)**2
+ (self.z-other.z)**2)

class Particle(Point):
"""The derived class for a particle."""
def __init__(self, coords=[], domain=1.0, m=1.0):
Point.__init__(self, coords, domain)
self.m = m # Mass of the particle (charge for potential)
self.phi = 0.  # Gravitational potential, initialized for direct sum

def P2M(sources, center, p):
"""
Given sources and cell's information, return multipole moments.
This computes M_n^m = sum q_i r_i^n * conj(Y_n^m(theta_i, phi_i))
"""
xc, yc, zc = center.x, center.y, center.z

# Coordinates of sources relative to the center
x_rel = np.array([source.x for source in sources]) - xc
y_rel = np.array([source.y for source in sources]) - yc
z_rel = np.array([source.z for source in sources]) - zc
masses = np.array([source.m for source in sources])

r_rel = np.sqrt(x_rel**2 + y_rel**2 + z_rel**2)
phi_rel = np.arctan2(y_rel, x_rel)

# Handle sources exactly at the center (r_rel = 0)
# np.arccos(0/0) results in NaN, set theta to 0 for points at center
theta_rel = np.arccos(z_rel / np.where(r_rel == 0, 1e-10, r_rel)) # Avoid division by zero
theta_rel[r_rel == 0] = 0.0 # Standard convention for (0,0,0) in spherical coords

multipole_moments = np.zeros((p + 1, 2 * p + 1), complex)

for i in range(len(sources)): # Loop over each source particle
for n in range(p + 1):
for m in range(-n, n + 1):
# q_i * r_i^n * conj(Y_n^m(theta_i, phi_i))
multipole_moments[n, m + n] += masses[i] * (r_rel[i]**n) * np.conj(sph_harm(m, n, phi_rel[i], theta_rel[i]))
return multipole_moments

def eval_potential(targets, multipole_moments, center, p):
"""
Given targets list, multipole moments and expansion center, return
the array of target's potentials.
This evaluates Phi(r') = sum (4pi/(2n+1)) * M_n^m / (r')^(n+1) * Y_n^m(theta', phi')
"""
potential_at_targets = np.zeros(len(targets), dtype=complex)

for j, target in enumerate(targets):
# Coordinates of target relative to the expansion center
dx = target.x - center.x
dy = target.y - center.y
dz = target.z - center.z

r_target = np.sqrt(dx**2 + dy**2 + dz**2)

# This expansion is valid only if r_target > r_max_sources_in_cell.
# For a target exactly at the center (r_target=0), the expansion is singular.
# In a full FMM, P2P (direct) interaction handles targets close to sources/expansion center.
if r_target < 1e-10: # If target is too close to center, expansion is invalid
potential_at_targets[j] = np.nan # Indicate invalid result for this case
continue

phi_target = np.arctan2(dy, dx)
theta_target = np.arccos(dz / r_target)

current_potential_sum = 0.0j # Initialize sum for this specific target
for n in range(p + 1):
for m in range(-n, n + 1):
# Y_n^m(theta', phi')
harmonic_val = sph_harm(m, n, phi_target, theta_target)

# Term for current (n, m)
term = (4 * np.pi / (2 * n + 1)) * (1 / (r_target**(n + 1))) * harmonic_val * multipole_moments[n, m + n]
current_potential_sum += term
potential_at_targets[j] = current_potential_sum

return potential_at_targets.real # Potentials are real

def direct_sum(sources, targets):
"""
Calculate the gravitational potential (target.phi) at each target
particle using direct summation method.
"""
phi_values = np.zeros(len(targets))
for j, target in enumerate(targets):
current_phi = 0.0
for source in sources:
r = target.distance(source)
if r < 1e-10: # Handle near-singularities / self-interaction
current_phi += np.inf # For gravitational potential, it goes to infinity
else:
current_phi += source.m / r
phi_values[j] = current_phi
return phi_values

# --- Test Setup for M2M ---
if __name__ == "__main__":
n_sources = 10  # number of particles in a child cell
p = 8         # order of expansion

# 1.  Define child cell sources and center
child_center_coords = [0.25, 0.25, 0.25]
child_center = Point(child_center_coords)

np.random.seed(42) # for reproducibility
# Sources distributed within a small radius around the child center
# Ensure they are within a sphere that defines the child cell
radius_child = 0.1
child_source_coords = np.random.uniform(-radius_child, radius_child,(n_sources, 3)) + np.array(child_center_coords)
child_sources = [Particle(coord.tolist(), m=1.0) for coord in child_source_coords]

# 2. Compute child's multipole moments (P2M)
child_multipole = P2M(child_sources, child_center, p)
print(f"Child multipole moments computed for {n_sources} sources.")

# 3. Define parent cell center
# Parent center should be "outside" the child cell, making d non-zero

parent_center_coords = [0.5, 0.5, 0.5] # E.g., a larger cell's center

parent_center = Point(parent_center_coords)
print(f"Parent center: {parent_center_coords}")
print(f"Child center: {child_center_coords}")

# 4. Perform M2M translation
parent_multipole_direct = P2M(child_sources, parent_center, p)

parent_multipole_M2M = M2M_theorem53(child_multipole, child_center, parent_center, p)

print(f"M2M translation performed.")

# 5. Define a distant target to test potentials
# The target must be sufficiently far from *both* child and parent centers
# for the multipole expansions to be accurate.
target_coords = [[2.0, 2.0, 2.0]] # Very distant target
targets = [Particle(coord) for coord in target_coords]
print(f"Target for evaluation: {target_coords[0]}")

# 6. Evaluate potential using child_multipole (centered at child_center)
# The target must be outside the smallest sphere enclosing the sources centered at child_center
phi_from_child_m = eval_potential(targets, child_multipole, child_center, p)
print(f"Potential from child multipole: {phi_from_child_m[0]:.12e}")

# 7. Evaluate potential using parent_multipole (centered at parent_center)
# The target must be outside the smallest sphere enclosing the sources centered at parent_center
phi_from_parent_m = eval_potential(targets, parent_multipole_M2M, parent_center, p)
print(f"Potential from parent multipole: {phi_from_parent_m[0]:.12e}")

# 8. Direct sum for comparison (ground truth)
# This involves directly summing contributions from all original child sources to the target
phi_direct = direct_sum(child_sources, targets)
print(f"Potential from direct sum: {phi_direct[0]:.12e}")

# 9. Check agreement
error_child_vs_direct = np.abs(phi_from_child_m[0] - phi_direct[0])
error_parent_vs_direct = np.abs(phi_from_parent_m[0] - phi_direct[0])
error_m2m_consistency = np.abs(phi_from_parent_m[0] - phi_from_child_m[0])

print(f"\nErrors (absolute difference):")
print(f"  Child multipole vs Direct: {error_child_vs_direct:.2e}")
print(f"  Parent multipole vs Direct: {error_parent_vs_direct:.2e}")
print(f"  M2M consistency (Parent vs Child potential): {error_m2m_consistency:.2e}")

if error_m2m_consistency < 1e-9: # A small tolerance for floating point comparisons
print("\nM2M translation appears to be working correctly: "
"Potential from parent's moments matches child's moments.")
else:
print("\nWarning: M2M translation might have an issue. "
"Check the consistency error.")

Das Ergebnis, das ich erhalten habe, ist:

Code: Select all

Child multipole moments computed for 10 sources.
Parent center: [0.5, 0.5, 0.5]
Child center: [0.25, 0.25, 0.25]
M2M translation performed.
Target for evaluation: [2.0, 2.0, 2.0]
Potential from child multipole: 3.275435467795e+00
Potential from parent multipole: 1.123673185858e+00
Potential from direct sum: 3.275435467795e+00

Errors (absolute difference):
Child multipole vs Direct: 8.26e-14
Parent multipole vs Direct: 2.15e+00
M2M consistency (Parent vs Child potential): 2.15e+00

Warning: M2M translation might have an issue. Check the consistency error.