Exercise_02/molecular_dynamics.py

@@ -88,7 +88,7 @@ def energy(position, box_size):
     # enforce expected shape
     # (`scipy.optimize.minimize` drops shape information)
     if len(position.shape) == 1:
-        position = position.reshape((-1, 3))
+        position = position.reshape((position.shape[0] // 3, 3))
     # compute all pairwise position differences (all!)
     diff = position[:, jnp.newaxis, :] - position
     # extract only one of two distance combinations of non-equal particles
@@ -122,16 +122,16 @@ def kinetic(velocity):
     return (mass / 2.0) * (velocity**2).sum()
 
 @jit
-def step(position, velocity, acceleration, box_size, delta_t):
+def step(position, velocity, acceleration, box_size, step_size):
     """
-    Performs a single Newton time step with `delta_t` given system state
+    Performs a single Newton time step with `step_size` given system state
     through the current particle `position`, `velocity` and `acceleration`.
     """
     # update position with a second order Taylor expantion
-    position += delta_t * velocity + (0.5 * delta_t**2) * acceleration
+    position += step_size * velocity + (0.5 * step_size**2) * acceleration
     # compute new particle acceleration through Newton’s second law of motion
     acceleration_next = force(position, box_size) / mass
     # update velocity with a finite mean approximation
-    velocity += (0.5 * delta_t) * (acceleration + acceleration_next)
+    velocity += (0.5 * step_size) * (acceleration + acceleration_next)
     # updated state
     return position, velocity, acceleration_next

Exercise_02/task02.py

@@ -59,8 +59,7 @@ from molecular_dynamics import dump, energy, force, mass
 position = np.random.uniform(0.0, config.box_size, (config.nr_particles, 3))
 
 # Sample particle velocities
-K_b = scipy.constants.Boltzmann / scipy.constants.eV # [eV / K]
+sd = np.sqrt(scipy.constants.Boltzmann * config.temperature / mass)
-sd = np.sqrt(K_b * config.temperature / mass)
 velocity = np.random.normal(0.0, sd, (config.nr_particles, 3))
 # center velocities
 velocity -= velocity.mean(axis = 0)

Exercise_03/NSSC_1.m

@@ -1,42 +1,28 @@
-%% Euler forward Dirichlet
+%% eulero forward DIrichlet
-L = 1;  % domain size (in the lecture notes this is denote h)
+L=pi;
-T = 2;  % time limit (max time)
+T=2;
-f = @(x,t) 0*x.*t;    % rhs of the more general equation `u_t - d u_xx = f`
+f=@(x,t) 0*x.*t;
-c1 = @(t) 1+0*t;      % _right_ boundary condition
+c1=@(t) 1+0*t;
-c2 = @(t) 0*t;        % _left_ boundary condition
+c2=@(t) 0*t;
-u0 = @(x) 0*x;        % initial values
+u0=@(x) 0*x;
-D = 0.5;              % diffusion parameter `d` in `u_t - d u_xx = f`
+D=0.5;
-%uex = @(x,t) cos(x).*exp(t);
+%uex=@(x,t) cos(x).*exp(t);
 
-N = 10;               % nr. of _space_ discretization points
+N=10;
-K = 200;              % nr. of _time_ discretization points
+K=200;
-[x, t, u] = Dirichlet_EA(L, N, T, K, c1, c2, f, u0, D);
+[x,t,u]=Dirichlet_EA(L,N,T,K,c1,c2,f,u0,D);
-
-% Report stability condition `D Delta T / (Delta x)^2 > 0.5`
-Delta_T = T / K;
-Delta_x = L / N;
-d = D * Delta_T / Delta_x^2;
-fprintf("Stability Condition: 0.5 >= D * Delta_T / Delta_x^2 = %f\n", d)
-if d > 0.5
-    fprintf("-> NOT Stable\n")
-else
-    fprintf("-> Stable\n")
-end
 
 figure(1)
-for ii = 1:K+1 % iterates time
+for ii=1:K+1 
-    hold on
+    plot(x,u(:,ii)');
-    plot(x, u(:, ii)');
     xlim([0 L])
     pause(0.05);
-    hold off
 end
 
-% 3D plot of space solution over time
+space=linspace(0,L,101);
-space = linspace(0,L,101);
+time=linspace(0,T,201);
-time = linspace(0,T,201);
+[xx,yy]=meshgrid(time,space);
-[xx,yy] = meshgrid(time,space);
+%exsol=uex(yy,xx);
-%exsol = uex(yy,xx);
 figure(2)
 mesh(t,x,u)
 %figure(2)

Exercise_03/task01_1-2.py

@@ -1,63 +0,0 @@
-# Task 1.1, 1.2
-import numpy as np
-from typing import Callable
-from matplotlib import pyplot as plt
-
-# Config
-D = 1e-6            # diffusion coefficient
-h = 1               # space domain (max x size)
-T = 2e6             # solution end time
-nx = 50             # nr of space discretization points
-nt = 20000          # nr of time discretization points
-
-# derived constants
-dx = h / (nx - 1)   # space step size
-dt = T / (nt - 1)   # time step size
-d = dt * D / dx**2  # stability/stepsize coefficient
-
-# report stability
-if d > 0.5:
-    print("NOT Stable")
-else:
-    print("Stable")
-
-# explicit scheme integration for `u_t - D u_xx = 0` with boundary conditions
-# enforced by `set_bounds` and initial conditions `initial`.
-def integrate(*, name: str, initial: np.array, set_bounds: Callable[[np.array], None]) -> None:
-    C = initial
-    # Setup boundary conditions
-    set_bounds(C)
-
-    i = 0       # index for plot generation
-    plt.figure(figsize = (8, 6), dpi = 100)
-    for t in range(nt):
-        # every 400'th time step save a plot
-        if t % (nt // 400) == 0:
-            plt.clf()
-            plt.plot(np.linspace(0, h, nx), C)
-            plt.xlim([0, h])
-            plt.ylim([0, 1.2])
-            plt.savefig(f"plots/{name}_{i:0>5}.png")
-            i += 1
-        # update solution using the explicit schema
-        C[1:-1] += d * (C[2:] - 2 * C[1:-1] + C[:-2])
-        # update right Neumann BC
-        set_bounds(C)
-
-# Subtask 1 boundary conditions (Dirichlet and Neumann)
-def bounds_1(C):
-    C[0] = 1
-    C[-1] = C[-2]
-
-# Subtask 2 boundary conditions (two Dirichlet)
-def bounds_2(C):
-    C[0] = 1
-    C[-1] = 0
-
-# run simulations
-integrate(name = 'task01_1', initial = np.zeros(nx), set_bounds = bounds_1)
-integrate(name = 'task01_2', initial = np.zeros(nx), set_bounds = bounds_2)
-
-# to convert generated image sequence to video use:
-# $> ffmpeg -r 60 -i plots/task01_1_%05d.png -pix_fmt yuv420p video_1_1.mp4
-# $> ffmpeg -r 60 -i plots/task01_2_%05d.png -pix_fmt yuv420p video_1_2.mp4

Exercise_03/task01_3.py

@@ -1,50 +0,0 @@
-# Task 1.3
-import numpy as np
-from typing import Callable
-from matplotlib import pyplot as plt
-
-# Config
-D = 1e-6            # diffusion coefficient
-h = 1               # space domain (max x size)
-T = 2e6             # solution end time
-nx = 50             # nr of space discretization points
-nt = 20000          # nr of time discretization points
-
-# derived constants
-dx = h / (nx - 1)   # space step size
-dt = T / (nt - 1)   # time step size
-d = dt * D / dx**2  # stability/stepsize coefficient
-
-# Setup implicit scheme equation matrix
-T = (1 + 2 * d) * np.eye(nx) - d * np.eye(nx, k = 1) - d * np.eye(nx, k = -1)
-# fix boundary condition equations
-T[0, 0] = 1         # Left Dirichlet BC
-T[0, 1] = 0
-T[-1, -2] = 1       # Right Neumann BC
-T[-1, -1] = 0
-
-# Set initial solution
-C = np.zeros(nx)
-C[0] = 1
-C[-1] = C[-2] # (0 = 0)
-
-i = 0       # index for plot generation
-plt.figure(figsize = (8, 6), dpi = 100)
-for t in range(nt):
-    # every 400'th time step save a plot
-    if t % (nt // 400) == 0:
-        plt.clf()
-        plt.plot(np.linspace(0, h, nx), C)
-        plt.xlim([0, h])
-        plt.ylim([0, 1.2])
-        plt.savefig(f"plots/task01_3_{i:0>5}.png")
-        i += 1
-    # update solution using the explicit schema
-    C = np.linalg.solve(T, C)
-    # fix BC conditions (theoretically, they are set by the update but for
-    # stability reasons (numerical) we enforce the correct values)
-    C[0] = 1
-    C[-1] = C[-2]
-
-# to convert generated image sequence to video use:
-# $> ffmpeg -r 60 -i plots/task01_3_%05d.png -pix_fmt yuv420p video_1_3.mp4