Recap
import numpy as np
import matplotlib.pyplot as pltRecap en uitbreiding¶
Je hebt al gezien hoe je een botsend deeltjesmodel maakt. Je startte toen met 1 of 2 deeltjes en kan dit uitbreiden naar ~100 deeltjes. Daarbij is het logisch om gebruik te maken van classes.
Voor het maken van een simulatie van gasmoleculen in een gesloten volume zal je logischerwijs de volgende stappen doorlopen:
bepalen van initiële condities
definiëren van de deeltjes
definiëren van volume met randvoorwaarden
simulaties bestaande uit:
bepaling van positie van deeltjes
controleren op onderlinge botsingen
controleren op botsingen met wanden
aanpassen van eigenschappen zoals totaal volume of temperatuur
opslaan van gegevens
De initiële condities¶
In onze simulatie gaan we uit van een aantal aannames. Bij deze aannames wil je dat de simulatie een voorspellende waarde heeft voor een reëel systeem, maar ook dat de benodigde rekenkracht voor het uitvoeren van de simulatie op te brengen is. We kiezen daarom voor de volgende condities:
Onze model heeft een beperkt en constant aantal deeltjes.
We maken gebruik van een 2D simulatie.
De deeltjes voelen geen onderlinge kracht en ondergaan alleen elastische botsingen.
(tijdens simuleren) De tijdstap is voldoende klein. (d.w.z. de afgelegde weg is klein in vergelijking met de diameter van de deeltjes, zodat we geen botsingen missen)
Hieronder geven we een aantal constanten voor onze code. Om de simulatie straks te kunnen interpreteren, vergelijken we deze eerst met een realistische situatie:
De simulatie die we gaan maken bevat maar twee dimensies. We kunnen het oppervlak per deeltje in deze simulatie daarom kiezen als het zojuist berekende .
Voor de gemiddelde snelheid, , van een molecuul geldt dat
Hierbij is de massa van het deeltje, het aantal dimensies, de constante van Boltzmann en de temperatuur.
Onze simulatie modelleert dus een klein systeem zonder rekening te houden met de kracht tussen de luchtmoleculen. Je zult zien dat we toch al behoorlijk wat van de processor vragen om deze berekeningen te doen. Om een meer realistische modelsimulatie te maken is daarom een serieuze uitdaging, die veel wiskunde en programmeerervaring vraagt. Op dit ogenblik wordt dit bijvoorbeeld gevraagd bij een van de eerste opdrachten van het vak ‘Computational Physics’ van de master Applied Physics.
# Ruimte voor uitwerking
BOX_SIZE_0 = 0 # Hoogte en breedte startvolume (in vul je gekozen eenheid in)
N = 40 # Aantal deeltjes
V_0 = 0 # Startsnelheid van deeltjes (vul je gekozen eenheid in)
RADIUS = 0 # Straal van moleculen (vul je gekozen eenheid in)
DT = 0 # Tijdstap om geen botsing te missen (vul je gekozen eenheid in)
### begin-solution
BOX_SIZE_0 = 10 # Hoogte en breedte startvolume (nm)
N = 40 # Aantal deeltjes
V_0 = 1 # Startsnelheid van deeltjes (in nm per 2.5 ps)
RADIUS = 0.3 # Straal van moleculen (nm)
DT = 0.1 * RADIUS / V_0 # Tijdstap om geen botsing te missen (* 2.5 ps)
### end-solutionParticle class¶
Ons ParticleClass definieert van een deeltje de massa , snelheid , positie , en straal . De positie moet per tijdstap bepaald worden.
Nieuw aan ParticleClass zijn twee properties, nl. momentum en kin_energy. Een property is eigenlijk een variabele die samenhangt met de waarde van andere variabelen binnen de klasse. Om de interne consistentie te bewaren moet je de juiste waarde van deze variabele dus telkens afleiden en dat doe je met een functie. Zo geeft kin_energy de kinetische energie terug van het deeltje door die af te leiden van de variabelen m en v.
Buiten ParticleClass voegen we twee functies toe, die je ook al hebt gezien. De eerste heet collide_detection en detecteert of twee deeltjes onderling overlappen en dus botsen. De tweede heet particle_collision en berekent de nieuwe snelheden van de twee deeltjes ten gevolge van hun botsing.
# Maken van de class
class ParticleClass:
def __init__(self, m, v, r, R):
""" maakt een deeltje (constructor) """
self.m = m
self.v = np.array(v, dtype=float)
self.r = np.array(r, dtype=float)
self.R = R
def update_position(self):
""" verandert positie voor één tijdstap """
self.r += self.v * DT
@property
def momentum(self):
return self.m * self.v
@property
def kin_energy(self):
return 1/2 * self.m * np.dot(self.v, self.v)
def collide_detection(p1: ParticleClass, p2: ParticleClass) -> bool:
""" Geeft TRUE als de deeltjes overlappen """
""" Geeft TRUE als de deeltjes overlappen """
dx = p1.r[0]-p2.r[0]
dy = p1.r[1]-p2.r[1]
rr = p1.R + p2.R
return dx**2+dy**2 < rr**2
#Note: np.linalg.norm(p1.r - p2.r) < (p1.R + p2.R) had ook gekund, maar is veel langzamer. Controleer zelf!
def particle_collision(p1: ParticleClass, p2: ParticleClass):
""" past snelheden aan uitgaande van overlap """
m1, m2 = p1.m, p2.m
delta_r = p1.r - p2.r
delta_v = p1.v - p2.v
dot_product = np.dot(delta_r, delta_v)
# Als deeltjes van elkaar weg bewegen dan geen botsing
if dot_product >= 0: # '='-teken voorkomt ook problemen als delta_r == \vec{0}
return
distance_squared = np.dot(delta_r, delta_r)
# Botsing oplossen volgens elastische botsing in 2D
p1.v -= 2 * m2 / (m1 + m2) * dot_product / distance_squared * delta_r
p2.v += 2 * m1 / (m1 + m2) * dot_product / distance_squared * delta_rVolume en randvoorwaarden¶
Voor het volume maken we gebruik van het oppervlak van een pyplot, net zoals we dat in de vorige werkbladen gedaan hebben. Daar is op dit ogenblik niets anders voor nodig dan een functie die de deeltjes laat botsen met de wanden.
Je hebt hiervoor al een voorbeeld van gezien, maar hieronder gebruiken we een andere vorm. Als je code ontwikkelt is het zaak om regelmatig en structureel te controleren of je code klopt. Het maakt daarbij niet uit of de code door jezelf is gemaakt, door een ander (zoals de docent) is aangeboden of van een generatieve AI komt. Hoe langer je code toevoegt zonder een controle uit te voeren hoe groter de kans is dat je fout op fout stapelt en niet meer kan traceren wat er precies mis is.
We maken gebruik van een truc die in veel simulaties wordt gebruikt: we maken gebruik van periodieke randvoorwaarden. In dit geval wil dit zeggen dat een deeltje dat botst met de wand aan de ene kant van het volume, aan de andere kant verschijnt met exact dezelfde snelheid. Op deze manier verandert de snelheid van een deeltje niet bij het botsen met de wand en blijven alle behoudswetten gelden. Die behoudswetten kunnen we straks goed controleren.
def box_collision(particle: ParticleClass):
''' botsing met wanden geeft periodieke randvoorwaarden '''
for i in range(2):
if particle.r[i] > BOX_SIZE_0 / 2:
particle.r[i] -= BOX_SIZE_0
elif particle.r[i] < -BOX_SIZE_0 / 2:
particle.r[i] += BOX_SIZE_0Maken van de simulatie¶
Nu dat we alle benodigde functies hebben, kunnen we ze in de juiste structuur zetten om de simulatie uit te voeren.
Bepaling positie deeltjes¶
Eerst moeten we alle deeltjes aanmaken en de startwaardes van hun variabelen invoeren:
def create_particles(particles):
''' Leegmaken en opnieuw aanmaken van deeltjes '''
particles.clear
for i in range(N):
vx = np.random.uniform(-V_0, V_0)
vy = np.random.choice([-1, 1]) * np.sqrt(V_0**2 - vx**2)
pos = np.random.uniform(-BOX_SIZE_0/2 + RADIUS, BOX_SIZE_0/2 - RADIUS, 2)
particles.append(ParticleClass(m=1.0, v=[vx, vy], r=pos, R=RADIUS))Controleren op onderlinge botsingen¶
Om de botsingen van alle deeltjes goed te verwerken moeten we controleren welke paren deeltjes met elkaar moeten botsen en deze botsingen één keer uitvoeren.
def handle_collisions(particles):
""" alle onderlinge botsingen afhandelen voor deeltjes in lijst """
num_particles = len(particles)
for i in range(num_particles):
for j in range(i+1, num_particles):
if collide_detection(particles[i], particles[j]):
particle_collision(particles[i], particles[j])Controleren op botsingen met wanden¶
Voor het bepalen of deeltjes botsen met de wanden maken we gebruik van de functie die we net hebben aangepast.
def handle_walls(particles):
''' botsing met wanden controleren voor alle deeltjes in volume '''
for p in particles:
box_collision(p)Integreren in tijdstap¶
Nu moeten we alle functies op de juiste manier combineren om de simulatie een tijdstap te laten maken:
def take_time_step(particles):
for p in particles:
p.update_position()
handle_collisions(particles)
handle_walls(particles)
Uitvoeren simulatie en tonen resultaat¶
Zoals aangegeven, kunnen we een animatie maken van de positie en snelheid als functie van de tijd, maar we kunnen ook het eindresultaat tonen en interpreteren:
particles = []
create_particles(particles)
for i in range(100):
take_time_step(particles)
plt.figure()
plt.xlabel('x')
plt.ylabel('y')
plt.gca().set_aspect('equal')
plt.xlim(-BOX_SIZE_0/2, BOX_SIZE_0/2)
plt.ylim(-BOX_SIZE_0/2, BOX_SIZE_0/2)
for p in particles:
plt.plot(p.r[0], p.r[1], 'k.', ms=25)
plt.arrow(p.r[0], p.r[1], p.v[0], p.v[1],
head_width=0.05, head_length=0.1, color='red')
plt.show()Toevoegen animatie¶
Zoals gezegd: het kost veel meer rekenkracht om een animatie te maken over de tijd. Soms kan het wel helpen om te zien of de simulatie aan alle verwachtingen voldoet.
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
fig, ax = plt.subplots()
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_xlim(-BOX_SIZE_0/2, BOX_SIZE_0/2)
ax.set_ylim(-BOX_SIZE_0/2, BOX_SIZE_0/2)
ax.set_aspect('equal')
dot, = ax.plot([], [], 'ro', ms=13);
def init():
dot.set_data([], [])
return dot,
def update(frame):
take_time_step(particles)
dot.set_data([p.r[0] for p in particles], [p.r[1] for p in particles])
return dot,
ani = FuncAnimation(fig, update, frames=100, init_func=init, blit=True, interval=50)
HTML(ani.to_jshtml())Controle code¶
Om te controleren of onze code werkt passen we nu de code aan.
# plotten van kinetische energie als functie van tijd
particles = []
energy_tot = []
create_particles(particles)
for i in range(100):
take_time_step(particles)
# voeg de som van de kinetische energie in deze stap toe aan energy_tot
### begin-solution
energy_tot.append(sum(p.kin_energy for p in particles))
### end-solution
# plot de figuur van de kinetische energie als functie van de tijd.
### begin-solution
plt.figure()
plt.xlabel('$t$')
plt.ylabel('$E_kin$')
plt.xlim(-0, 100)
plt.ylim(N/2 - 1e-9, N/2 + 1e-9)
plt.plot(np.linspace(1, 100, 100), energy_tot, '-b')
plt.show()
### end-solution### begin-solution
particles = []
px_tot = []
create_particles(particles)
for i in range(100):
take_time_step(particles)
px_tot.append(sum(p.momentum[0] for p in particles))
plt.figure()
plt.xlabel('$t$')
plt.ylabel('$p_x$')
plt.xlim(-0, 100)
plt.ylim(min(px_tot),max(px_tot))
plt.plot(np.linspace(1, 100, 100), px_tot, '-b')
plt.show()
### end-solutionVoor simulaties waarbij de kans op botsingen heel groot is (richting vloeistof), is er een handige methode om de snelheid van je simulatie op te voeren. In onderstaande code wordt daar een voorbeeld van gegeven.
def handle_collisions(particles):
ignore_list = []
for i, p1 in enumerate(particles):
if p1 in ignore_list:
continue
for j, p2 in enumerate(particles):
if p1 is p2:
continue
if collide_detection(particles[i], particles[j]):
particle_collision(particles[i], particles[j])
ignore_list.append(p2)In de volgende notebooks maken we veelal gebruik van numpy functies. In sommige gevallen is dat sneller, maar lang niet altijd. Het is good-practice om na te gaan wat sneller is. Aan de andere kant moet de tekst ook leesbaar blijven...
Als onderdeel om excellent te scoren verwachten we dat je een aantal checks tussendoor doet om code sneller te maken. In onderstaande testcode, gemaakt door Josh Sleijfer wordt een van die functies getest. Wat blijkt? Voor “kleine” vectoren is het sneller om geen gebruik te maken van Numpy maar voor grotere vectoren wel!
import numpy as np
import time
L = 2 # lengte van de vector
repeats = 1_000_000 # aantal herhalingen
rnd_vec1 = np.random.rand(L, repeats)
# test numpy dot product
t0 = time.perf_counter()
for vec1 in rnd_vec1.T:
dot = np.dot(vec1, vec1)
t1 = time.perf_counter()
print(f"Numpy: Time taken for {repeats} dot products of length {L}: {t1 - t0:.4f} seconds")
# test no numpy dot product
t0 = time.perf_counter()
for vec1 in rnd_vec1.T:
dot = vec1[0] * vec1[0] + vec1[1] * vec1[1] # + vec1[2] * vec1[2] + vec1[3] * vec1[3] # don't forget to adjust!
t1 = time.perf_counter()
print(f"No Numpy: Time taken for {repeats} dot products of length {L}: {t1 - t0:.4f} seconds")
# test fully vectorized
t0 = time.perf_counter()
dot = np.sum(rnd_vec1[0] * rnd_vec1[0] + rnd_vec1[1] * rnd_vec1[1])
t1 = time.perf_counter()
print(f"Fully vectorized: Time taken for {repeats} dot products of length {L}: {t1 - t0:.4f} seconds")