Project 6 by Group 2

projects/project-6-group-2.ipynb

@@ -0,0 +1,169 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib notebook\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "from matplotlib import animation\n",
+    "from IPython.display import HTML"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def grad_LJ(r, eps, sig):\n",
+    "    n = len(r)\n",
+    "    grad = np.zeros_like(r, dtype=\"float64\")\n",
+    "    for i in range(0,n):\n",
+    "        for j in range(0, n):\n",
+    "            if i != j:\n",
+    "                dr = r[i] - r[j]\n",
+    "                length = np.linalg.norm(dr)\n",
+    "                grad[i] += eps*((6*sig**6)/length**8 - (12*sig**12)/length**12)*dr\n",
+    "    return grad "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def grad_harmonic(r, k):\n",
+    "    return 2*k*r"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def potential_gradient(r, eps=1, sig=1, k=1):\n",
+    "    return (grad_LJ(r, eps, sig) + grad_harmonic(r, k))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def vv(potential_gradient, r_init, v_init, dt=1/100, t_ges=100, eps=1., sig=1., k=1.):\n",
+    "    \n",
+    "    size = int(t_ges/dt)\n",
+    "    n = len(r_init)\n",
+    "    \n",
+    "    r_matrix, v_matrix = np.zeros((size, n, 2)), np.zeros((size, n, 2))\n",
+    "    r_matrix[0], v_matrix[0] = r_init, v_init\n",
+    "    \n",
+    "    for i in range(1,size):\n",
+    "        \n",
+    "        r = r_matrix[i-1]\n",
+    "        v = v_matrix[i-1]\n",
+    "        gradient = potential_gradient(r, eps, sig, k)\n",
+    "        \n",
+    "        r_new = r + dt * v - 1/2 * dt**2 * gradient\n",
+    "        gradient_new = potential_gradient(r_new, eps, sig, k) \n",
+    "        v_new = v - 1/2 * dt * (gradient + gradient_new)\n",
+    "        \n",
+    "        r_matrix[i], v_matrix[i] = r_new, v_new\n",
+    "\n",
+    "        \n",
+    "    return r_matrix, v_matrix"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "r_init = np.random.uniform(low=-5, high=5,size=(5,2))\n",
+    "v_init = np.zeros_like(r_init)\n",
+    "r_matrix, v_matrix = vv(potential_gradient, r_init, v_init, k=1/100)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots(figsize=(5, 5))\n",
+    "scat = ax.scatter(r_matrix[0,:,0], r_matrix[0,:,1], c=np.arange(len(r_init)))\n",
+    "qax = ax.quiver(r_matrix[0,:,0], r_matrix[0,:,1], v_matrix[1,:,0], v_matrix[1,:,1],np.arange(len(r_init)),scale=20, width=0.005)\n",
+    "ax.set_xlim(np.max([-10,np.min(r_matrix[:,:,0])]),np.min([10,np.max(r_matrix[:,:,0])]))\n",
+    "ax.set_ylim(np.max([-10,np.min(r_matrix[:,:,1])]),np.min([10,np.max(r_matrix[:,:,1])]))\n",
+    "\n",
+    "def animate(i):\n",
+    "    index = 4*i\n",
+    "    data = r_matrix[index]\n",
+    "    scat.set_offsets(data)\n",
+    "    qax.set_UVC(v_matrix[index,:,0],v_matrix[index,:,1])\n",
+    "    qax.set_offsets(data)\n",
+    "\n",
+    "anim = animation.FuncAnimation(fig, animate, interval=20, blit=True, repeat=False)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def gd(potential_gradient, r_init, gamma=1/100, size=10000, eps=1, sig=1, k=1):\n",
+    "    n = len(r_init)\n",
+    "    r_ = r_init\n",
+    "    for i in range(1,size):\n",
+    "        r = r_ -gamma*potential_gradient(r_, eps, sig, k)\n",
+    "        r_ = r[::]\n",
+    "    return r  "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "r_init = np.random.uniform(low=-1, high=1,size=(3,2))\n",
+    "r = gd(potential_gradient,r_init, k=1/10)\n",
+    "plt.scatter(r[:,0],r[:,1])\n",
+    "plt.xlim(-2,2)\n",
+    "plt.ylim(-2,2)\n",
+    "plt.show()\n",
+    "    "
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

projects/project_4/poisson_2.py

@@ -0,0 +1,48 @@
+import numpy as np
+
+def create_laplacian_2d(nx, ny, lx, ly, pbc=True):
+    """ Computes discrete Laplacian for a 2d
+        charge density matrix, ordered row-wise
+        Args:
+            nx(int  >= 2): number of grid points along x axis
+            ny(int  >= 2): number of grid points along y axis
+            lx(float > 0): length of grid along x axis
+            ly(float > 0): length of grid along y axis
+            pbc(boolean): periodic boundry conditions
+        output:
+            Laplacian as nx * ny by nx * ny np.array
+    """
+    if type(nx) != int or type(ny) != int:
+        raise TypeError('We need an integer')
+    if type(lx) != int and type(lx) != float:
+        raise TypeError('We need a number')
+    if type(ly) != int and type(ly) != float:
+        raise TypeError('We need a number')
+    if nx < 2 or ny < 2:
+        raise ValueError('We need at least two grid points')
+    if lx <= 0 or ly <= 0:
+        raise ValueError('We need a positive length')
+    if type(pbc) != bool:
+        raise TypeError('We need a boolean as pbc')
+
+    hx = (nx / lx) ** 2
+    hy = (ny / ly) ** 2
+    a1 =  np.diag((-2 * hx - 2 * hy) * np.ones(nx * ny))
+    diag1 = hx * np.ones(nx * ny - 1)
+    diag1[nx-1::nx] = 0
+    a2 = np.diag(diag1 , 1)
+    a3 = np.diag(diag1 , -1)
+    a4 = np.diag(hy * np.ones(nx * ny - nx), nx)
+    a5 = np.diag(hy * np.ones(nx * ny - nx), -nx)
+    laplacian = a1 + a2 + a3 + a4 + a5
+
+    if pbc:
+        a6 = np.diag(hy * np.ones(nx), nx * ny - nx)
+        a7 = np.diag(hy * np.ones(nx), -nx * ny + nx)
+        diag2 = hx * np.zeros(nx * ny - nx + 1)
+        diag2[::nx] = hx
+        a8 = np.diag(diag2, nx - 1)
+        a9 = np.diag(diag2, -nx + 1)
+        laplacian += a6 + a7 + a8 + a9
+
+    return laplacian

projects/project_5/jacobi.py

@@ -0,0 +1,58 @@
+import numpy as np
+
+def jacobi(charge_dist, boxsize, maxiter=1000, tol=10**(-5), epsilon0=1.0):
+
+    """An implementation of the Jacobi-method for solving the periodic discretized Poisson-equation in n dimensions.
+    Returns the potential as a n-dimensional array.
+
+    Arguments:
+        charge_dist: n-dimensional array consisting of entries of type float that define the charge distribution.
+        boxsize: tuple of length n, which containts the lengths of the box-axes.
+        maxiter: int value that sets the maximum number of iterations being performed by the algorithm.
+        tol: float that sets the break-off-condition for the algorithm.
+        epsilon0: float, scaling factor of the charge distribution.
+    """
+
+    if type(boxsize) in (int, float):
+        boxsize = (boxsize, )
+
+    h_array = np.asarray(boxsize)
+    rho = np.asarray(charge_dist)
+    pot = np.zeros(rho.shape)
+    dim = len(h_array)
+
+    if type(maxiter) is not int:
+        raise TypeError("Value maxiter has to be of type int.")
+    if type(tol) not in (int, float):
+        raise TypeError("Value tol has to be a real numeric value.")
+    if type(epsilon0) not in (int, float):
+        raise TypeError("Value epsilon0 has to be a real numeric value.")
+    if rho.dtype not in (int, float):
+        raise TypeError("Entries of charge_dist have to be of type int or float.")
+    if h_array.dtype not in (int, float):
+        raise TypeError("Entries of boxsize have to be either of type int or float.")
+    if h_array.ndim != 1:
+        raise TypeError("Boxsize needs to be a one-dimensional array.")
+    if dim != rho.ndim:
+        raise ValueError("Boxsize needs to have as many values as charge_dist has axes.")
+    if any(h_array==0.):
+        raise ValueError("Boxsize values all have to be nonzero.") 
+    if (1 or 0) in rho.shape:
+        raise ValueError("All axes of charge_dist need to have more than one element.")
+    if not np.isclose(np.mean(rho), 0, atol=1e-16):
+        raise ValueError("The mean of charge_dist has to be zero. Consider subtracting the mean and try again.")
+
+    c = 2*np.sum(1/(h_array)**2)
+    dif = np.inf
+    counter = 0  
+
+    while tol < dif and counter < maxiter:
+        pot_next = np.zeros(rho.shape)
+        for i, h in enumerate(h_array[::-1]):
+            pot_next += (np.roll(pot, 1, axis=i)+np.roll(pot, -1, axis=i))/h**2  
+        pot_next = (1/c)*(pot_next + rho/epsilon0)
+        dif = np.linalg.norm(pot_next - pot)
+        pot = pot_next
+        counter+= 1
+
+    return pot

projects/project_5/test/test_jacobi.py

@@ -0,0 +1,60 @@
+import pytest
+import numpy as np
+from itertools import product
+from project_5.jacobi import jacobi
+from project_4.poisson_2 import create_laplacian_2d
+
+@pytest.mark.parametrize('rho, boxsize, exception', [
+    ([[1,2,-3],[4,5,'yes']], (1,3), TypeError),
+    ([1,-1], None, TypeError),
+    ([1,2,-3], 'hello', TypeError),
+    ([1,2,-3], [[2, 1]], TypeError),
+    ([1,-1], 0, ValueError),
+    ([1,2,-3], [0,0.1], ValueError),
+    ([[1,-1],[-1,1]], (0.1,0.1,0.1), ValueError),
+    (-1, (2), ValueError),
+    ([[]], (0.1,0.1), ValueError),
+    ([[1,2,-3]],(0.1,0.1), ValueError),
+    ([1, 1], (0.1), ValueError)])
+def test_jacobi_exceptions(rho, boxsize, exception):
+    with pytest.raises(exception):
+        jacobi(rho, boxsize)
+
+@pytest.mark.parametrize('nx, ny, lx, ly', [
+    (nx, ny, lx, ly) for nx, ny, lx, ly in product(
+        [5, 10, 20],
+        [5, 10, 20],
+        [1.0, 3.0],
+        [1.0, 3.0])])
+def test_consisteny(nx, ny, lx, ly):
+    x = np.linspace(0, lx, nx)
+    y = np.linspace(0, ly, ny)
+    hx, hy = lx/nx, ly/ny
+    boxsize = (hx, hy)
+
+    rho = np.random.normal(size=(ny,nx))
+    rho -= np.mean(rho)
+    laplacian = create_laplacian_2d(nx, ny, lx, ly)
+    pot = jacobi(rho, boxsize, maxiter=5000)
+
+    assert pot.ndim == 2, \
+    f'Potentials have wrong dimensions: {pot.ndim}'
+
+    assert pot.shape[0] == ny, \
+    f'Potentials have wrong first shape: {pot.shape[0]}'
+
+    assert pot.shape[1] == nx, \
+    f'Potentials have wrong second shape: {pot.shape[1]}'
+
+    laplacian = create_laplacian_2d(nx, ny, lx, ly)
+    pot_laplacian = np.linalg.solve(laplacian, -rho.reshape(nx*ny))
+    pot_laplacian = pot_laplacian.reshape((ny,nx))
+    pot_laplacian = pot_laplacian + (pot[0, 0] - pot_laplacian[0, 0]) #Equivalent potentials can differ by a constant.
+
+    np.testing.assert_allclose(pot_laplacian, pot, rtol=1e-2, atol=1e-2)
+
+
+
+
+
+