Skip to content
Projects

Fields and a simulation

Download notebookView on GitHubRun when the site was built.

The field of 20,000 random charges by Barnes-Hut tree walk, compared with the exact direct sum; the opening angle theta trades accuracy for speed.

import numpy as np
import tree_d_3d
from tree_d_3d import FlatQuadtree, Quad, direct_field
rng = np.random.default_rng(0)
n = 20_000
pos = rng.uniform(-1.0, 1.0, (2, n)) # component-major: pos[0] is x
charge = rng.choice([-1.0, 1.0], n)
others = np.arange(n) # each particle excludes its own charge
reference = direct_field(pos, pos, charge, softening=1e-3, exclude=others)[1]
tree = FlatQuadtree(pos, charge, Quad(0.0, 0.0, 1.0))
for theta in (1.0, 0.5, 0.25):
field = tree.evaluate(pos, exclude=others, theta=theta, softening=1e-3)[1]
error = np.abs(field - reference).max() / np.abs(reference).max()
print(f"theta = {theta}: max. relative error {error:.1e}")
theta = 1.0: max. relative error 8.6e-03
theta = 0.5: max. relative error 1.5e-04
theta = 0.25: max. relative error 6.5e-06

For line charges (kernel="log", the 2D Poisson kernel) the fast multipole method computes the fields of all particles at once in O(N).

lines = FlatQuadtree(pos, charge, Quad(0.0, 0.0, 1.0), kernel="log", order=16)
reference = direct_field(pos, pos, charge, exclude=others, kernel="log")[1]
field = lines.fmm(theta=0.5)[1]
print(
f"FMM: max. relative error {np.abs(field - reference).max() / np.abs(reference).max():.1e}"
)
FMM: max. relative error 3.1e-10

An ElectrostaticSimulation advances ChargedBody particles; a cloud of like charges expands, converting potential into kinetic energy.

bodies = [tree_d_3d.ChargedBody(x, y) for x, y in rng.normal(0.0, 0.1, (500, 2))]
sim = tree_d_3d.ElectrostaticSimulation(
bodies, Quad(0.0, 0.0, 1.0), k_e=1.0, softening=1e-3, dt=1e-4
)
energy = [sim.get_total_energy()]
for _ in range(200):
sim.step()
energy.append(sim.get_total_energy())
print(f"relative energy change {abs(energy[-1] - energy[0]) / abs(energy[0]):.1e}")
relative energy change 1.3e-03