Python Fast Multiplopel реализация с использованием сферических гармоник ⇐ Python

[Anonymous]

Я пытаюсь реализовать метод Fast Multipole после короткого курса
http://math.nyu.edu/faculty/greengar/sh ... se_fmm.pdf.
Теперь предположим, что поставленная задача состоит в том, чтобы вычислить взаимодействие кулонов в широкой системе части. Если точечная плата за прочность $ Q $ расположена по адресу $ p_0 = \ left (x_0, y_0, z_0 \ right) $, то потенциал и электростатическое поле из -за этой зарядки в отдельной точке $ p = (x, y, z) $ даются < /p>
$$ \ phi = \ frac {1 {\ tild {r}} {r}}} {r}}}} {r}} < /p> $ $ $ \ tilde {r} $ обозначает расстояние между точками $ p_0 $ и $ p $. < /p>
Теорема 5.2 (Multiople расширение) Предположим, что $ k $ заряды сильных сторон $ \ Left {q_i, i = 1, \ ldots, k \ right} $, расположенные в оставшихся точках \ \ beta_i \ right), i = \ right. $ 1, \ ldots, k} $, с $ \ left | \ rho_i \ right | a $, потенциал $ \ phi (p) $ определяется уравнением (5-15): < /p>
$$ \ phi (p) = \ sum_ {n = 0 {\ {\ infty} \ sum_ {n = -n^{n = -n {n = -n {n = {n = {n = {n = {n \ \ \ frac {m_n^m} {r^{n+1}} \ cdot y_n^m (\ theta, \ phi) $$ < /p>
где < /p>
$$ m_n^m = \ sum_ {i = 1}^k a \ cdot \ rho_i^n \ Y_n^{-m} \ left (\ alpha_i, \ beta_i \ right) $$ < /p>
Теорема 5.3 (перевод мультипольного расширения) Предположим, что $ l $ заряды за сильные стороны $ Q_1, Q_2, \ cdots, Q_L $ расположены внутри $ D $ D $ a at q = wh h h h h ah ah ah h h \ beta) $, и что для баллов $ p = (r, \ theta, \ phi) $ werld $ d $, потенциал из-за этих сборов определяется мультиполиным расширением < /p>
$$ \ phi (p) = \ sum_ {n = 0}^{\ infty} \ sum_ {m = -n}^{rac {\ {{\ {{\ {{O_N^{{O_N^{{O_ {{{O_N^{raC} n+1}} \ cdot y_n^m \ left (\ theta^{\ prime}, \ phi^{\ prime} \ right) $$ < /p>
где $ p-q = \ left (r^{\ prime}, \ theta^{\ prime}, \ phi^{\ prime} \ right) $. Тогда для любой точки $ p = (r, \ theta, \ phi) $ за пределами сферы $ d_1 $ of radius $ (a+\ rho) $, < /p>
$$ \ phi (p) = \ sum_ {j = 0}^{\ infty} \ sum_ {k = -j}^j \ frac {m_j^k} {1 k} {1 k} {1 k^k}}} {1 k} \ cdot y_j^k (\ theta, \ phi) $$
и
$$ m_j^k = \ sum_ {n = 0}^j \ sum_ {m = -n}^n \ frac {o_ {j-n}^{k-m} \ cdot i^{| k |-| m |-| k-m |} {k-m} \ cdot i^{| A_ {j-n}^{k-m} \ cdot \ rho^n \ cdot y_n^{-m} (\ alpha, \ beta)} {a_j^k} $$
с $ a_n^m $ defined,
$ a_n^m $ /> $$ a_n^m = \ frac {(-1)^n} {\ sqrt {(n-m)! (n+m)!}} $$ < /p>
и < /p>
$$ o_n^m = \ sum_ {i = 1}^n q_i \ cdot r_i r_i r_i^n r_i Y_n^{-m} \ left (\ theta_i, \ phi_i \ right) $$ < /p>
$ o_n^m $ = Old Multipolole (Multiople расширение уже рассчитано). Его центр находится в детском центре. Мультиполь вычислен с использованием родительского центра, но с использованием уже вычисленного $ o_n^m $ рассматривает детский центр. /> определение сферической гармоники из Sympy < /h1>
$$ y_n^m (\ theta, \ varphi): = \ sqrt {\ frac {(2 n+1) (n-m)!} {4 \ pi (n+m)!} \ exp (i m \ varphi) \ mathrm {n+m)} \ \ exp (i m \ varphi) \ \ pi (n+m)}}}}}}}}}}}} \ exp (i m \ varphi) {4 \ pi (n+m)! (\ theta)) $$ < /p>
определение сферической гармоники из короткого курса < /h1>
$$ y_n^m (\ theta, \ phi) = \ sqrt {\ frac {(n- | m |)!} {n+| m |)! \ theta) e^{i m \ phi} $$ < /p>
Обратите внимание, что в этом случае нет $ \ frac {(2 n+1)} {4 \ pi} $ в корне. Поскольку это различие, этот фактор должен появляться в реализации уравнения (5-15). Это в функции eval_potential в строке:
term = (4 * np.pi / (2 * n + 1)) * (1 / (r_target**(n + 1))) * harmonic_val * multipole_moments[n, m + n]
Моя реализация в порядке для мультипольного расширения, заданного уравнением (5-15) как потенциал из Multipolle: 3,275435467795e+00 очень близок к потенциалу от прямой суммы: 3.275435467795e+00. Однако перевод мультиполя, заданный уравнением (5-20), не работает должным образом, поскольку потенциал от родительского мультиполя: 1.123673185858E+00 сильно отличается от потенциала от прямой суммы. Ошибка находится в вычислении parent_multipole_direct. Значения переменной parent_multipole_direct должны быть равны значениям parent_multipole_m2m, но они отличаются. Я проверил весь код, но не могу найти ошибку. Может кто -нибудь мне помочь? Вот мой код: < /p>
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 * (r_rel**n) * np.conj(sph_harm(m, n, phi_rel, theta_rel))
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.")

< /code>
Результат, который я получил: < /p>
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.


Подробнее здесь: https://stackoverflow.com/questions/797 ... -harmonics