Created
November 29, 2023 21:00
-
-
Save ninux/4e44e845caed96e4a7fe2343c3e348c4 to your computer and use it in GitHub Desktop.
frist rudimentary FDTD GUI
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import threading | |
import time | |
import tkinter as tk | |
from tkinter import messagebox | |
import matplotlib.pyplot as plt | |
import numpy as np | |
from matplotlib.figure import Figure | |
from matplotlib.backends.backend_tkagg import (FigureCanvasTkAgg, | |
NavigationToolbar2Tk) | |
from matplotlib import animation | |
from matplotlib import style | |
class fdtd_core: | |
def __init__(self): | |
# create event to signal computation done | |
self.fdtd_new_data = threading.Event() | |
# general definitions and constants | |
self.h = 6.626e-34 # [J*s] Plank's constant | |
self.hbar = self.h / (2 * np.pi) # [J*s] reduced Plank's constant | |
self.m0 = 9.1e-31 # [kg] free space mass of electron (rest mass) | |
self.meff_Si = 1.08 # [-] effective mass factor of Si | |
self.meff_Ge = 0.067 # [-] effective mass factor of Ge | |
self.meff_GaAs = 0.55 # [-] effective mass factor of GaAs | |
self.meff = self.meff_Si # [m] effective mass | |
self.melec = self.meff * self.m0 # [?] mass of an electron | |
self.ecoul = 1.6e-19 # [C] charge of an electron | |
self.epsz = 8.85e-9 # [A*s/V/m] dielectric of free space | |
self.eV2J = 1.6e-19 # [-] energy conversion factor (eV to J) | |
self.J2eV = 1 / self.eV2J # [-] energy conversion factor (J to eV) | |
# energies | |
self.PE = float(0) | |
self.KE = float(0) | |
# simulation parameters | |
self.NN = 400 # [-] number of points in the problem space | |
self.del_x = 0.1e-9 # [?] cell size | |
self.dt = 2e-17 # [s] time steps | |
self.DX = self.del_x * 1e9 # [_] cell size in nm | |
self.XX = np.arange(0, self.DX * self.NN, self.DX) # [-] length in nm for plotting | |
self.ra = (0.5 * self.hbar / self.melec) * (self.dt / np.power(self.del_x, 2)) # must be < 0.15 for stability | |
# define stability threshold | |
self.ra_stability_threshold = 0.15 | |
if (self.ra < self.ra_stability_threshold): | |
# stable | |
pass | |
else: | |
# unstable | |
print("WARNING: Simulation setting results in unstable results") | |
self.V = np.zeros(self.NN, dtype=float) | |
# initialize sine wave in a gaussian envelope | |
self.pulse_lambda = 40 # pulse wavelength | |
self.pulse_sigma = 50 # pulse width | |
self.nc = 100 # starting position | |
self.ptot = 0.0 # total energy | |
self.prl = np.zeros(self.NN, dtype=float) # real part of the state variable | |
self.pim = np.zeros(self.NN, dtype=float) # imaginary part of the state variable | |
for n in range(1, self.NN - 1): | |
self.prl[n] = (np.exp(-1 * np.power((n-self.nc) / self.pulse_sigma, 2)) | |
* np.cos(2 * np.pi * (n-self.nc) / self.pulse_lambda)) | |
self.pim[n] = (np.exp(-1 * np.power((n-self.nc) / self.pulse_sigma, 2)) | |
* np.sin(2 * np.pi * (n-self.nc) / self.pulse_lambda)) | |
self.ptot = self.ptot + np.power(self.prl[n], 2) + np.power(self.pim[n], 2) | |
self.pnorm = np.sqrt(self.ptot) | |
# normalize data | |
for n in range(1, self.NN): | |
self.prl[n] = self.prl[n] / self.pnorm | |
self.pim[n] = self.pim[n] / self.pnorm | |
self.ptot = self.ptot + np.power(self.prl[n], 2) + np.power(self.pim[n], 2) | |
self.XX_copy = self.XX | |
self.prl_copy = self.prl | |
self.pim_copy = self.pim | |
self.T = 0 # absolute time in steps | |
self.n_step = 100 # steps per loop | |
self.count = 0 # number of loops | |
self.nos = 100 # number of timesteps between plots | |
self.nop = 3 # number of plots to be made | |
def fdtd_calculation(self): | |
for m in range(0, self.n_step - 1): | |
self.T += 1 | |
for n in range(1, self.NN - 1): | |
self.prl[n] = (self.prl[n] | |
- self.ra * (self.pim[n-1] - (2*self.pim[n]) + self.pim[n+1]) | |
+ (self.dt / self.hbar) * self.V[n] * self.pim[n]) | |
for n in range(1, self.NN - 1): | |
self.pim[n] = (self.pim[n] | |
+ self.ra * (self.prl[n-1] - (2*self.prl[n]) + self.prl[n+1]) | |
- (self.dt/self.hbar) * self.V[n] * self.prl[n]) | |
# copy data | |
self.XX_copy = self.XX | |
self.prl_copy = self.prl | |
self.pim_copy = self.pim | |
# set event so signal computation is completed, new data available | |
self.fdtd_new_data.set() | |
# calculate energies | |
ptot = float(0) | |
#PE = float(0) | |
psi = (self.prl*float(0) + 1j*self.pim*float(0)) | |
for n in range(0, self.NN): | |
psi[n] = self.prl[n] + 1j*self.pim[n] | |
self.PE = np.real(self.PE + (psi[n] * np.conjugate(psi[n]) * self.V[n])) | |
self.PE = self.PE * self.J2eV | |
ke = float(0) + 1j*float(0) | |
for n in range(1, self.NN-1): | |
lap_p = psi[n+1] - 2*psi[n] + psi[n-1] | |
ke = ke + lap_p*np.conjugate(psi[n]) | |
self.KE = -self.J2eV*( np.power(self.hbar/self.del_x, 2) / (2*self.melec)) * np.real(ke) | |
def get_location(self): | |
return self.XX_copy | |
def get_real(self): | |
return self.prl_copy | |
def get_imag(self): | |
return self.pim_copy | |
def fdtd_plot(self): | |
plt.plot(self.XX, self.prl, 'k-', self.XX, self.pim, 'k--') | |
plt.show() | |
def check_normalization(self): | |
self.norm_limit = (1e-3) * np.array([-1, 1]) + 1 | |
if (self.ptot > self.norm_limit[0] and self.ptot < self.norm_limit[1]): | |
# normalization ok | |
pass | |
else: | |
# normalization not ok | |
print("WARNING: Normalization for energy not ok") | |
def check_stability(self): | |
if (self.ra < self.ra_stability_threshold): | |
# stable | |
pass | |
else: | |
# unstable | |
print("WARNING: Simulation not stable, check parameters.") | |
def set_simulation_size(self, value): | |
self.NN = int(value) | |
def get_simulation_size(self): | |
return int(self.NN) | |
def set_steps(self, value): | |
self.n_step = int(value) | |
def get_steps(self): | |
return self.n_step | |
def get_potential_energy(self): | |
return self.PE | |
def get_kinetic_enregy(self): | |
return self.KE | |
def get_total_energy(self): | |
return (self.PE + self.KE) | |
def reset_simulation(self): | |
# initialize sine wave in a gaussian envelope | |
self.pulse_lambda = 40 # pulse wavelength | |
self.pulse_sigma = 50 # pulse width | |
self.nc = 100 # starting position | |
self.ptot = float(0.0) # total energy | |
self.prl = np.zeros(self.NN, dtype=float) # real part of the state variable | |
self.pim = np.zeros(self.NN, dtype=float) # imaginary part of the state variable | |
for n in range(1, self.NN - 1): | |
self.prl[n] = (np.exp(-1 * np.power((n - self.nc) / self.pulse_sigma, 2)) | |
* np.cos(2 * np.pi * (n - self.nc) / self.pulse_lambda)) | |
self.pim[n] = (np.exp(-1 * np.power((n - self.nc) / self.pulse_sigma, 2)) | |
* np.sin(2 * np.pi * (n - self.nc) / self.pulse_lambda)) | |
self.ptot = self.ptot + np.power(self.prl[n], 2) + np.power(self.pim[n], 2) | |
self.pnorm = np.sqrt(self.ptot) | |
self.ptot = float(0) | |
# normalize data | |
for n in range(0, self.NN): | |
self.prl[n] = self.prl[n] / self.pnorm | |
self.pim[n] = self.pim[n] / self.pnorm | |
self.ptot = self.ptot + np.power(self.prl[n], 2) + np.power(self.pim[n], 2) | |
print("ptot = " + str(self.ptot)) | |
class fdtdgui: | |
def __init__(self): | |
self.root = tk.Tk() | |
#style.use("dark_background") | |
# initialize FDTD simulation | |
self.fdtd = fdtd_core() | |
# initial geometry | |
self.root.geometry("900x500") | |
# parameters | |
self.plot_auto_number_value = 1 | |
# menubar | |
self.menubar = tk.Menu(self.root) | |
# close menu | |
self.closemenu = tk.Menu(self.menubar, tearoff=0) | |
self.closemenu.add_command(label="Close", command=self.on_closing) | |
self.menubar.add_cascade(menu=self.closemenu, label="File") | |
# about menu | |
self.aboutmenu = tk.Menu(self.menubar, tearoff=0) | |
self.aboutmenu.add_command(label="About", command=self.show_about) | |
self.menubar.add_cascade(menu=self.aboutmenu, label="About") | |
self.root.config(menu=self.menubar) | |
# plot figure | |
self.frame = tk.Frame(self.root) | |
self.fig, self.ax = plt.subplots() | |
self.canvas = FigureCanvasTkAgg(self.fig, master=self.frame) | |
self.frame.grid(row=0, column=0, rowspan=10, columnspan=1, sticky="nsew") | |
self.canvas.get_tk_widget().pack() | |
# program name label | |
self.name_label_frame = tk.Frame(self.root) | |
self.name_label = tk.Label(self.name_label_frame, | |
text="FDTD Simulation", | |
font=('Arial', 12, 'bold')) | |
self.name_label_frame.grid(row=0, column=1, columnspan=2, sticky="nsew") | |
self.name_label.pack(fill=tk.BOTH) | |
# plot single button | |
self.plot_single_button_frame = tk.Frame(self.root) | |
self.plot_single_button = tk.Button(self.plot_single_button_frame, | |
text="Plot Single", | |
font=('Arial', 10, 'bold'), | |
command=self.fdtd_plot_single) | |
self.plot_single_button_frame.grid(row=1, column=1, columnspan=2, sticky="nsew") | |
self.plot_single_button.pack(fill=tk.BOTH) | |
# plot auto button | |
self.plot_auto_button_frame = tk.Frame(self.root) | |
self.plot_auto_button = tk.Button(self.plot_auto_button_frame, | |
text="Plot Auto", | |
font=('Arial', 10, 'bold'), | |
command=self.fdtd_plot_auto) | |
self.plot_auto_button_frame.grid(row=2, column=1, columnspan=2, sticky="nsew") | |
self.plot_auto_button.pack(fill=tk.BOTH) | |
# plot reset button | |
self.plot_reset_button_frame = tk.Frame(self.root) | |
self.plot_reset_button = tk.Button(self.plot_reset_button_frame, | |
text="Reset Simulation", | |
font=('Arial', 10, 'bold'), | |
command=self.reset_simulation) | |
self.plot_reset_button_frame.grid(row=3, column=1, columnspan=2, sticky="nsew") | |
self.plot_reset_button.pack(fill=tk.BOTH) | |
# simulation parameters label | |
self.label_simulation_parameters_frame = tk.Frame(self.root) | |
self.label_simulation_parameters = tk.Label(self.label_simulation_parameters_frame, | |
text="Simulation Parameters", | |
font=('Arial', 10, 'bold')) | |
self.label_simulation_parameters_frame.grid(row=4, column=1, sticky="nsw") | |
self.label_simulation_parameters.pack(fill=tk.BOTH) | |
# steps label | |
self.steps_label_frame = tk.Frame(self.root) | |
self.steps_label = tk.Label(self.steps_label_frame, | |
text="Steps per plot (-):") | |
self.steps_label_frame.grid(row=5, column=1, sticky="nsw") | |
self.steps_label.pack(fill=tk.BOTH) | |
# steps slider | |
self.steps_slider_frame = tk.Frame(self.root) | |
self.steps_slider = tk.Scale(self.steps_slider_frame, | |
from_=1, to=1000, orient=tk.HORIZONTAL, | |
command=self.set_steps) | |
self.steps_slider.set(self.get_steps()) | |
self.steps_slider_frame.grid(row=5, column=2, sticky="nsew") | |
self.steps_slider.pack(fill=tk.BOTH) | |
# plot auto number label | |
self.plot_auto_number_label_frame = tk.Frame(self.root) | |
self.plot_auto_number_label = tk.Label(self.plot_auto_number_label_frame, text="Auto Plots (-):") | |
self.plot_auto_number_label_frame.grid(row=6, column=1, sticky="nsw") | |
self.plot_auto_number_label.pack(fill=tk.BOTH) | |
# plot auto number slider | |
self.plot_auto_number_slider_frame = tk.Frame(self.root) | |
self.plot_auto_number_slider = tk.Scale(self.plot_auto_number_slider_frame, | |
from_=1, to=100, orient=tk.HORIZONTAL, | |
command=self.set_plot_auto_number) | |
self.plot_auto_number_slider_frame.grid(row=6, column=2, sticky="nsew") | |
self.plot_auto_number_slider.pack(fill=tk.BOTH) | |
def on_closing(self): | |
if messagebox.askyesno(title="Quit?", message="Do you really want to quit?"): | |
self.root.destroy() | |
def show_about(self): | |
messagebox.showinfo(title="About", message="FDTD Simulation, Version 0.0") | |
def fdtd_plot_single(self): | |
calc = threading.Thread(target=self.fdtd.fdtd_calculation) | |
calc.start() | |
calc.join() | |
x = self.fdtd.get_location() | |
prl = self.fdtd.get_real() | |
pim = self.fdtd.get_imag() | |
ke = self.fdtd.get_kinetic_enregy() | |
pe = self.fdtd.get_potential_energy() | |
self.ax.clear() | |
self.ax.plot(x, prl, label="real") | |
self.ax.plot(x, pim, label="imag") | |
self.ax.grid(visible=True) | |
plt.ylim(-0.15, 0.15) | |
plt.legend(loc="upper right") | |
plt.xlabel("x (nm)") | |
plt.ylabel(r"$\Psi(x)$") | |
plt.title("KE = " + f"{ke:.3}" + " eV, PE = " + f"{pe:.3}" + " eV") | |
self.canvas.draw() | |
self.canvas.flush_events() | |
def fdtd_plot_auto(self): | |
for m in range(0, self.plot_auto_number_value): | |
self.fdtd_plot_single() | |
def set_steps(self, value): | |
self.fdtd.set_steps(int(value)) | |
def get_steps(self): | |
return self.fdtd.get_steps() | |
def set_plot_auto_number(self, value): | |
self.plot_auto_number_value = int(value) | |
def get_plot_auto_number(self): | |
return self.plot_auto_number_value | |
def reset_simulation(self): | |
self.fdtd.reset_simulation() | |
self.fdtd_plot_single() | |
app = fdtdgui() | |
tk.mainloop() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment