1#!/usr/bin/env python3
2"""
31D Vlasov-Poisson Solver
4
5This script implements a 1D electrostatic Vlasov-Poisson solver using
6operator splitting (Strang splitting) with cubic spline interpolation.
7Demonstrates plasma oscillations.
8
9Author: Plasma Physics Examples
10License: MIT
11"""
12
13import numpy as np
14import matplotlib.pyplot as plt
15from scipy.constants import e, m_e, epsilon_0, k
16from scipy.interpolate import CubicSpline
17from scipy.fft import fft, ifft, fftfreq
18
19
20class VlasovPoisson1D:
21 """
22 1D-1V Vlasov-Poisson solver using operator splitting.
23
24 ∂f/∂t + v∂f/∂x + (q/m)E∂f/∂v = 0
25 ∂E/∂x = ρ/ε₀ = (q/ε₀)(∫f dv - n₀)
26 """
27
28 def __init__(self, Nx, Nv, Lx, v_max, n0, T, q=-e, m=m_e):
29 """
30 Initialize solver.
31
32 Parameters:
33 -----------
34 Nx, Nv : int
35 Grid points in x and v
36 Lx : float
37 System length [m]
38 v_max : float
39 Maximum velocity [m/s]
40 n0 : float
41 Background density [m^-3]
42 T : float
43 Temperature [K]
44 q, m : float
45 Charge and mass
46 """
47 self.Nx = Nx
48 self.Nv = Nv
49 self.Lx = Lx
50 self.v_max = v_max
51 self.n0 = n0
52 self.T = T
53 self.q = q
54 self.m = m
55
56 # Spatial grid
57 self.x = np.linspace(0, Lx, Nx, endpoint=False)
58 self.dx = Lx / Nx
59
60 # Velocity grid
61 self.v = np.linspace(-v_max, v_max, Nv)
62 self.dv = 2 * v_max / Nv
63
64 # Distribution function f(x, v)
65 self.f = np.zeros((Nx, Nv))
66
67 # Electric field
68 self.E = np.zeros(Nx)
69
70 # Wavenumbers for FFT
71 self.kx = 2 * np.pi * fftfreq(Nx, Lx / Nx)
72
73 def initialize_maxwellian(self, perturbation=0.0, k_pert=1):
74 """
75 Initialize with perturbed Maxwellian.
76
77 f(x, v, t=0) = n₀(1 + α cos(kx)) * f_M(v)
78
79 Parameters:
80 -----------
81 perturbation : float
82 Amplitude of density perturbation
83 k_pert : int
84 Mode number of perturbation
85 """
86 v_th = np.sqrt(2 * k * self.T / self.m)
87
88 # Maxwellian in velocity
89 f_M = (self.n0 / (np.sqrt(2 * np.pi) * v_th) *
90 np.exp(-self.v**2 / (2 * v_th**2)))
91
92 # Spatial perturbation
93 k = 2 * np.pi * k_pert / self.Lx
94 n_pert = 1 + perturbation * np.cos(k * self.x)
95
96 # f(x, v)
97 for i in range(self.Nx):
98 self.f[i, :] = n_pert[i] * f_M
99
100 def compute_density(self):
101 """Compute density n(x) = ∫ f(x,v) dv."""
102 return np.trapz(self.f, self.v, axis=1)
103
104 def compute_electric_field(self):
105 """
106 Solve Poisson equation for E field using FFT.
107
108 ∂E/∂x = ρ/ε₀ = (q/ε₀)(n - n₀)
109 """
110 n = self.compute_density()
111 rho = self.q * (n - self.n0)
112
113 # FFT of charge density
114 rho_k = fft(rho)
115
116 # Solve in Fourier space: ik E_k = rho_k / ε₀
117 # E_k = -i * rho_k / (k * ε₀)
118 E_k = np.zeros_like(rho_k, dtype=complex)
119 E_k[1:] = -1j * rho_k[1:] / (self.kx[1:] * epsilon_0)
120 E_k[0] = 0 # Zero mode (charge neutrality)
121
122 # Inverse FFT
123 self.E = np.real(ifft(E_k))
124
125 def advect_x(self, dt):
126 """
127 Advection in x: ∂f/∂t + v∂f/∂x = 0
128
129 Use cubic spline interpolation for accuracy.
130 """
131 f_new = np.zeros_like(self.f)
132
133 for j in range(self.Nv):
134 # For each velocity, advect in x
135 # x_new = x - v * dt (backwards)
136 x_old = (self.x - self.v[j] * dt) % self.Lx # Periodic BC
137
138 # Interpolate
139 cs = CubicSpline(self.x, self.f[:, j], bc_type='periodic')
140 f_new[:, j] = cs(x_old)
141
142 self.f = f_new
143
144 def advect_v(self, dt):
145 """
146 Advection in v: ∂f/∂t + (q/m)E∂f/∂v = 0
147
148 Use cubic spline interpolation.
149 """
150 f_new = np.zeros_like(self.f)
151 a = self.q * self.E / self.m # Acceleration
152
153 for i in range(self.Nx):
154 # For each position, advect in v
155 # v_new = v - a * dt (backwards)
156 v_old = self.v - a[i] * dt
157
158 # Handle boundaries: extrapolate or zero
159 mask = (v_old >= -self.v_max) & (v_old <= self.v_max)
160
161 if np.any(mask):
162 cs = CubicSpline(self.v, self.f[i, :], bc_type='natural')
163 f_new[i, mask] = cs(v_old[mask])
164 f_new[i, mask] = np.maximum(f_new[i, mask], 0) # Non-negative
165
166 self.f = f_new
167
168 def step_strang(self, dt):
169 """
170 Strang splitting: A(dt/2) B(dt) A(dt/2)
171
172 A = x-advection, B = v-advection
173 """
174 # Half step in x
175 self.advect_x(dt / 2)
176
177 # Full step in v (requires E field)
178 self.compute_electric_field()
179 self.advect_v(dt)
180
181 # Half step in x
182 self.advect_x(dt / 2)
183
184 def run(self, t_end, dt):
185 """
186 Run simulation.
187
188 Parameters:
189 -----------
190 t_end : float
191 End time [s]
192 dt : float
193 Timestep [s]
194
195 Returns:
196 --------
197 history : dict
198 Time history of fields
199 """
200 n_steps = int(t_end / dt)
201
202 # Storage
203 t_history = np.zeros(n_steps + 1)
204 E_history = np.zeros((n_steps + 1, self.Nx))
205 n_history = np.zeros((n_steps + 1, self.Nx))
206 f_history = []
207
208 # Initial state
209 self.compute_electric_field()
210 t_history[0] = 0
211 E_history[0, :] = self.E
212 n_history[0, :] = self.compute_density()
213 f_history.append(self.f.copy())
214
215 # Time loop
216 for n in range(n_steps):
217 self.step_strang(dt)
218
219 t_history[n + 1] = (n + 1) * dt
220 E_history[n + 1, :] = self.E
221 n_history[n + 1, :] = self.compute_density()
222
223 # Store f occasionally
224 if n % max(1, n_steps // 20) == 0:
225 f_history.append(self.f.copy())
226
227 return {
228 't': t_history,
229 'E': E_history,
230 'n': n_history,
231 'f_snapshots': f_history,
232 'x': self.x,
233 'v': self.v
234 }
235
236
237def simulate_plasma_oscillation():
238 """Simulate plasma oscillation and verify frequency."""
239
240 # Parameters
241 n0 = 1e16 # m^-3
242 T = 1e4 # K (~ 1 eV)
243 v_th = np.sqrt(2 * k * T / m_e)
244
245 # Plasma frequency
246 omega_pe = np.sqrt(n0 * e**2 / (epsilon_0 * m_e))
247 f_pe = omega_pe / (2 * np.pi)
248 T_pe = 1 / f_pe
249
250 print(f"\nPlasma Parameters:")
251 print(f" Density: {n0:.2e} m^-3")
252 print(f" Temperature: {T:.2e} K ({T * k / e:.2f} eV)")
253 print(f" Thermal velocity: {v_th:.2e} m/s")
254 print(f" Plasma frequency: {f_pe:.2e} Hz")
255 print(f" Plasma period: {T_pe:.2e} s")
256
257 # Debye length
258 lambda_D = np.sqrt(epsilon_0 * k * T / (n0 * e**2))
259 print(f" Debye length: {lambda_D:.2e} m")
260
261 # Grid
262 Nx = 64
263 Nv = 128
264 Lx = 10 * lambda_D # System size
265 v_max = 6 * v_th
266
267 # Initialize solver
268 solver = VlasovPoisson1D(Nx, Nv, Lx, v_max, n0, T)
269 solver.initialize_maxwellian(perturbation=0.01, k_pert=1)
270
271 # Run simulation
272 t_end = 10 * T_pe
273 dt = T_pe / 50
274
275 print(f"\nSimulation:")
276 print(f" Grid: {Nx} × {Nv}")
277 print(f" System size: {Lx:.2e} m ({Lx/lambda_D:.1f} λ_D)")
278 print(f" Time: {t_end:.2e} s ({t_end/T_pe:.1f} T_pe)")
279 print(f" Timestep: {dt:.2e} s ({dt/T_pe:.3f} T_pe)")
280
281 history = solver.run(t_end, dt)
282
283 return history, omega_pe, lambda_D
284
285
286def plot_results(history, omega_pe, lambda_D):
287 """Plot simulation results."""
288
289 t = history['t']
290 E = history['E']
291 n = history['n']
292 x = history['x']
293 v = history['v']
294 f_snapshots = history['f_snapshots']
295
296 T_pe = 2 * np.pi / omega_pe
297
298 # Create figure
299 fig = plt.figure(figsize=(16, 12))
300
301 # Plot 1: Electric field vs time (at one location)
302 ax1 = plt.subplot(3, 3, 1)
303 i_mid = len(x) // 2
304 ax1.plot(t / T_pe, E[:, i_mid], 'b-', linewidth=2)
305 ax1.set_xlabel('Time [T_pe]', fontsize=11)
306 ax1.set_ylabel('E(x=L/2) [V/m]', fontsize=11)
307 ax1.set_title('Electric Field Oscillation', fontsize=12, fontweight='bold')
308 ax1.grid(True, alpha=0.3)
309
310 # Plot 2: E field energy vs time
311 ax2 = plt.subplot(3, 3, 2)
312 E_energy = np.sum(E**2, axis=1) * (x[1] - x[0])
313 ax2.semilogy(t / T_pe, E_energy, 'r-', linewidth=2)
314 ax2.set_xlabel('Time [T_pe]', fontsize=11)
315 ax2.set_ylabel('∫E² dx [a.u.]', fontsize=11)
316 ax2.set_title('Electric Field Energy', fontsize=12, fontweight='bold')
317 ax2.grid(True, alpha=0.3)
318
319 # Plot 3: Density perturbation
320 ax3 = plt.subplot(3, 3, 3)
321 n0 = np.mean(n[0, :])
322 dn = (n - n0) / n0
323 cs = ax3.contourf(x / lambda_D, t / T_pe, dn, levels=20, cmap='RdBu_r')
324 ax3.set_xlabel('x [λ_D]', fontsize=11)
325 ax3.set_ylabel('Time [T_pe]', fontsize=11)
326 ax3.set_title('Density Perturbation δn/n₀', fontsize=12, fontweight='bold')
327 plt.colorbar(cs, ax=ax3)
328
329 # Plot 4-9: Phase space snapshots
330 snapshot_indices = [0, len(f_snapshots)//4, len(f_snapshots)//2,
331 3*len(f_snapshots)//4, len(f_snapshots)-1]
332
333 for idx, snap_idx in enumerate(snapshot_indices[:6]):
334 ax = plt.subplot(3, 3, 4 + idx)
335 f = f_snapshots[snap_idx]
336
337 # Time of snapshot
338 t_snap = snap_idx * (t[-1] / (len(f_snapshots) - 1))
339
340 # Plot phase space
341 v_th = np.sqrt(2 * k * 1e4 / m_e)
342 cs = ax.contourf(x / lambda_D, v / v_th, f.T, levels=20, cmap='hot')
343 ax.set_xlabel('x [λ_D]', fontsize=10)
344 ax.set_ylabel('v [v_th]', fontsize=10)
345 ax.set_title(f't = {t_snap/T_pe:.2f} T_pe', fontsize=11, fontweight='bold')
346 ax.set_ylim([-3, 3])
347
348 plt.tight_layout()
349 plt.savefig('vlasov_plasma_oscillation.png', dpi=300, bbox_inches='tight')
350 plt.show()
351
352 # FFT analysis
353 fig2, ax = plt.subplots(figsize=(10, 6))
354
355 # FFT of E field at midpoint
356 E_mid = E[:, i_mid]
357 dt = t[1] - t[0]
358 freqs = fftfreq(len(t), dt)
359 E_fft = np.abs(fft(E_mid))
360
361 # Only positive frequencies
362 pos_freqs = freqs > 0
363 ax.semilogy(freqs[pos_freqs] / (omega_pe / (2*np.pi)), E_fft[pos_freqs],
364 'b-', linewidth=2)
365 ax.axvline(1.0, color='red', linestyle='--', linewidth=2,
366 label=f'ω_pe = {omega_pe/(2*np.pi):.2e} Hz')
367
368 ax.set_xlabel('Frequency [f_pe]', fontsize=12)
369 ax.set_ylabel('|FFT(E)| [a.u.]', fontsize=12)
370 ax.set_title('Frequency Spectrum of E Field', fontsize=13, fontweight='bold')
371 ax.legend(loc='best', fontsize=11)
372 ax.grid(True, alpha=0.3)
373 ax.set_xlim([0, 3])
374
375 plt.tight_layout()
376 plt.savefig('vlasov_frequency_spectrum.png', dpi=300, bbox_inches='tight')
377 plt.show()
378
379
380if __name__ == '__main__':
381 print("\n" + "="*80)
382 print("1D VLASOV-POISSON SOLVER: PLASMA OSCILLATIONS")
383 print("="*80)
384
385 history, omega_pe, lambda_D = simulate_plasma_oscillation()
386
387 print("\nGenerating plots...")
388 plot_results(history, omega_pe, lambda_D)
389
390 print("\nDone! Generated 2 figures:")
391 print(" - vlasov_plasma_oscillation.png")
392 print(" - vlasov_frequency_spectrum.png")