Dislocation_dynamics-v01.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dislocation dynamics (2D)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some python code to animate dislocation dynamics for teaching purposes. Needs to be tidied up, decribed and made more pythonic (and faster)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.animation as animation\n",
    "%matplotlib qt5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def initDD(n,D):\n",
    "#  we use the rand function, which creates a 1D array of random\n",
    "#  numbers between 0 and 1\n",
    "    x_array = np.random.random(n)*D\n",
    "    y_array = np.random.random(n)*D\n",
    "    b_array = np.ones(n)\n",
    "# assign b with equal numbers of +1 and -1\n",
    "    b_array[:n//2]=1\n",
    "    b_array[n//2:]=-1\n",
    "    alive=np.ones(n, dtype=bool)\n",
    "    return x_array,y_array,b_array,alive\n",
    "\n",
    "def sumDD(n,D,rc,x_pos,y_pos,b_list,alive):\n",
    "#  set force component to 0\n",
    "    dx_list=[]\n",
    "    dy_list=[]\n",
    "    for x,y in zip(x_pos,y_pos):\n",
    "        dx = x-x_pos\n",
    "        dy = y-y_pos\n",
    "        dx = dx - D*np.round(dx/D)\n",
    "        dy = dy - D*np.round(dy/D)\n",
    "        dx_list.append(dx)\n",
    "        dy_list.append(dy)\n",
    "\n",
    "        \n",
    "    fx_array=np.zeros(n)\n",
    "    for dis1 in range(n):\n",
    "        fx=np.zeros(n)\n",
    "        for dis2 in range(n):\n",
    "            dx=dx_list[dis1][dis2]\n",
    "            dy=dy_list[dis1][dis2]\n",
    "            b1=b_list[dis1]\n",
    "            b2=b_list[dis2]\n",
    "            dsq = dx**2 + dy**2\n",
    "            dist = np.sqrt(dsq)\n",
    "            if np.logical_and(alive[dis1], alive[dis2]):\n",
    "                if np.logical_and(dist<1.2,b1==-b2):\n",
    "                    alive[dis1]=False\n",
    "                    alive[dis2]=False\n",
    "            if np.logical_and(alive[dis1], alive[dis2]):\n",
    "                if (0.9 < dist <= rc):\n",
    "                    fx[dis1] = fx[dis1] + b1*b2*dx*(dx**2-dy**2)/dsq**2\n",
    "            fx_array=fx_array+fx\n",
    "\n",
    "# calculate maximum value of force\n",
    "    fmax=max(abs(fx_array))\n",
    "    return fx_array,fmax,alive\n",
    "\n",
    "#  input:  ndis = the number of dislocations\n",
    "#          nsteps = number of time steps\n",
    "#          dxmax = maximum move per time step\n",
    "#\n",
    "#  output: xi=initial x positions,x = final positions,y,b\n",
    "#          fx = final forces on each dislocation\n",
    "#          xdm = the maximum distance a dislocation has moved in the\n",
    "#          calculation\n",
    "def DD2D(ndis,nsteps,dxmax):\n",
    "    xdm_list=[]\n",
    "    dx_list=[]\n",
    "    x_pos_list=[]\n",
    "    y_pos_list=[]\n",
    "    s_list=[]\n",
    "    # set size of system (arbitrary)\n",
    "    D = 100\n",
    "    # set cutoff to be 1/2 the box length\n",
    "    rc = D/2\n",
    "    # initial positions\n",
    "    x_pos,y_pos,b_list,alive = initDD(ndis,D)\n",
    "    x_pos_list.append(x_pos.copy())\n",
    "    y_pos_list.append(y_pos.copy())\n",
    "    # store initial position\n",
    "    xi = x_pos.copy()\n",
    "    # start the time steps\n",
    "    for step in range(nsteps):\n",
    "        fx_array,fmax,alive = sumDD(ndis,D,rc,x_pos,y_pos,b_list,alive)\n",
    "        dt = dxmax/fmax\n",
    "        x_pos=x_pos+fx_array*dt\n",
    "        escr=np.where(x_pos>D)\n",
    "        escl=np.where(x_pos<0)\n",
    "        x_pos[escr]=x_pos[escr]-D\n",
    "        x_pos[escl]=x_pos[escl]+D\n",
    "        xd = x_pos - xi#change in position during run\n",
    "        xdm=np.max(abs(xd))# find maximum movement\n",
    "        xdm_list.append(xdm.copy())\n",
    "        dx_list.append(np.sum(np.abs(fx_array*dt)))\n",
    "        x_pos_list.append(x_pos.copy())\n",
    "        y_pos_list.append(y_pos.copy())\n",
    "        s_list.append(alive)\n",
    "        xd = xd - D*np.round(xd/D)# remove movement across periodic boundaries\n",
    "    return xi,x_pos_list,y_pos_list,b_list,xdm_list,dx_list,s_list\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sumDD_climb(n,D,rc,x_pos,y_pos,b_list,alive):\n",
    "#  set force component to 0\n",
    "    dx_list=[]\n",
    "    dy_list=[]\n",
    "    for x,y in zip(x_pos,y_pos):\n",
    "        dx = x-x_pos\n",
    "        dy = y-y_pos\n",
    "        dx = dx - D*np.round(dx/D)\n",
    "        dy = dy - D*np.round(dy/D)\n",
    "        dx_list.append(dx)\n",
    "        dy_list.append(dy)\n",
    "     \n",
    "    fx_array=np.zeros(n)\n",
    "    fy_array=np.zeros(n)\n",
    "    for dis1 in range(n):\n",
    "        fx=np.zeros(n)\n",
    "        fy=np.zeros(n)\n",
    "        for dis2 in range(n):\n",
    "            dx=dx_list[dis1][dis2]\n",
    "            dy=dy_list[dis1][dis2]\n",
    "            b1=b_list[dis1]\n",
    "            b2=b_list[dis2]\n",
    "            dsq = dx**2 + dy**2\n",
    "            dist = np.sqrt(dsq)\n",
    "            if np.logical_and(alive[dis1], alive[dis2]):\n",
    "                if np.logical_and(dist<1.5,b1==-b2):\n",
    "                    alive[dis1]=False\n",
    "                    alive[dis2]=False\n",
    "            if np.logical_and(alive[dis1], alive[dis2]): # check for anihilation\n",
    "                if (0.9 < dist <= rc):\n",
    "                    fx[dis1] = fx[dis1] + b1*b2*dx*(dx**2-dy**2)/dsq**2\n",
    "                    fy[dis1] = fy[dis1] - b1*b2*dy*(3*dx**2+dy**2)/dsq**2\n",
    "            fx_array=fx_array+fx\n",
    "            fy_array=fy_array+fy\n",
    "    # calculate maximum value of force\n",
    "    fmax=max([max(abs(fx_array)),max(abs(fy_array))])\n",
    "    return fx_array, fy_array ,fmax,alive\n",
    "\n",
    "#  input:  ndis = the number of dislocations\n",
    "#          nsteps = number of time steps\n",
    "#          dxmax = maximum move per time step\n",
    "#\n",
    "#  output: xi=initial x positions,x = final positions,y,b\n",
    "#          fx = final forces on each dislocation\n",
    "#          xdm = the maximum distance a dislocation has moved in the\n",
    "#          calculation\n",
    "def DD2D_climb(ndis,nsteps,dxmax):\n",
    "    xdm_list=[]\n",
    "    ydm_list=[]\n",
    "    dy_list=[]\n",
    "    dx_list=[]\n",
    "    x_pos_list=[]\n",
    "    y_pos_list=[]\n",
    "    s_list=[]\n",
    "    # set size of system (arbitrary)\n",
    "    D = 100\n",
    "    # set cutoff to be 1/2 the box length\n",
    "    rc = D/2\n",
    "    # initial positions\n",
    "    x_pos,y_pos,b_list,alive = initDD(ndis,D)\n",
    "    x_pos_list.append(x_pos.copy())\n",
    "    y_pos_list.append(y_pos.copy())\n",
    "    s_list.append(alive.copy())\n",
    "    # store initial position\n",
    "    xi = x_pos.copy()\n",
    "    yi = y_pos.copy()\n",
    "    # start the time steps\n",
    "    for step in range(nsteps):\n",
    "        fx_array,fy_array,fmax,alive = sumDD_climb(ndis,D,rc,x_pos,y_pos,b_list,alive)\n",
    "        dt = dxmax/fmax\n",
    "        x_pos=x_pos+fx_array*dt\n",
    "        y_pos=y_pos+fy_array*dt*0.1\n",
    "        escr=np.where(x_pos>D)\n",
    "        escl=np.where(x_pos<0)\n",
    "        x_pos[escr]=x_pos[escr]-D\n",
    "        x_pos[escl]=x_pos[escl]+D\n",
    "        escr=np.where(y_pos>D)\n",
    "        escl=np.where(y_pos<0)\n",
    "        y_pos[escr]=y_pos[escr]-D\n",
    "        y_pos[escl]=y_pos[escl]+D\n",
    "        xd = x_pos - xi#change in position during run\n",
    "        yd= x_pos - yi #where is yi?\n",
    "        #continue from here\n",
    "        xdm=np.max(abs(xd))# find maximum x movement\n",
    "        xdm_list.append(xdm.copy())\n",
    "        ydm=np.max(abs(yd))\n",
    "        ydm_list.append(ydm.copy())\n",
    "        dx_list.append(np.sum(np.abs(fx_array*dt)))\n",
    "        dy_list.append(np.sum(np.abs(fy_array*dt)))\n",
    "        x_pos_list.append(x_pos.copy())\n",
    "        y_pos_list.append(y_pos.copy())\n",
    "        s_list.append(alive.copy())\n",
    "        xd = xd - D*np.round(xd/D)# remove movement across periodic boundaries\n",
    "        yd = yd - D*np.round(yd/D)\n",
    "    return xi,x_pos_list,y_pos_list,b_list,xdm_list,dx_list,s_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_dis(x_pos,y_pos,n):\n",
    "    mbot = ((-7, 0), (0, 0), (0, -8), (0, 0), (7, 0))\n",
    "    mtop = ((-7, 0), (0, 0), (0, 8), (0, 0), (7, 0))\n",
    "    fig, ax = plt.subplots()\n",
    "    scat = ax.scatter(x_pos[:n//2],y_pos[:n//2],c='red',s=200,cmap='viridis',marker=mtop)\n",
    "    scat2 = ax.scatter(x_pos[n//2:],y_pos[n//2:],c='blue',s=200,cmap='viridis',marker=mbot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "mbot = ((-7, 0), (0, 0), (0, -8), (0, 0), (7, 0))\n",
    "mtop = ((-7, 0), (0, 0), (0, 8), (0, 0), (7, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test functions (no climb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initiate dislocation array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_pos_list=[]\n",
    "y_pos_list=[]\n",
    "for n in range(10):\n",
    "    D =100\n",
    "    rc=D/2\n",
    "    n_dis=20\n",
    "    x_pos,y_pos,b_list,alive = initDD(n_dis,D)\n",
    "    #plot_dis(x_pos,y_pos,n_dis)\n",
    "    x_pos_list.append(x_pos.copy())\n",
    "    y_pos_list.append(y_pos.copy())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create position lists with climb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_dis=20\n",
    "n_steps=400\n",
    "max_dx=0.3\n",
    "xi,x_pos_list_glide,y_pos_list_glide,b_list,xdm,f_res,s_list_glide = DD2D(n_dis,n_steps,max_dx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create positions with climb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n_dis=20\n",
    "n_steps=500\n",
    "max_dx=0.7\n",
    "xi,x_pos_list_climb,y_pos_list_climb,b_list,xdm,f_res,s_list_glide = DD2D_climb(n_dis,n_steps,max_dx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Animate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Choose glide only or climb:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_pos_list = x_pos_list_glide\n",
    "x_pos_list = x_pos_list_glide\n",
    "s_list = s_list_glide\n",
    "\n",
    "#x_pos_list = x_pos_list_climb\n",
    "#x_pos_list = x_pos_list_glide\n",
    "#s_list = s_list_glide"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "mbot = ((-7, 0), (0, 0), (0, -8), (0, 0), (7, 0))\n",
    "mtop = ((-7, 0), (0, 0), (0, 8), (0, 0), (7, 0))\n",
    "\n",
    "D=100\n",
    "n=n_dis\n",
    "fig = plt.figure(figsize=(10,10))\n",
    "ax=plt.axes(xlim=(0,D),ylim=(0,D))\n",
    "line_up, = ax.plot([],[],marker=mtop,ms=20,lw=0,mew=2)\n",
    "line_down, = ax.plot([],[],marker=mbot,ms=20,lw=0,mew=2)\n",
    "time_step =ax.text(0,101,'')\n",
    "\n",
    "\n",
    "def init():\n",
    "    line_up.set_data([],[]) \n",
    "    line_down.set_data([],[])\n",
    "    time_step.set_text('Time step = 0')             \n",
    "    return line_up, line_down, time_step\n",
    "    \n",
    "def animate(i):\n",
    "    x_up=x_pos_list[i][:n//2]\n",
    "    x_up=x_up[np.where(s_list[i][:n//2])]\n",
    "    y_up=y_pos_list[i][:n//2]\n",
    "    y_up=y_up[np.where(s_list[i][:n//2])]\n",
    "    x_down=x_pos_list[i][n//2:]\n",
    "    x_down=x_down[np.where(s_list[i][n//2:])]\n",
    "    y_down=y_pos_list[i][n//2:]\n",
    "    y_down=y_down[np.where(s_list[i][n//2:])]\n",
    "    line_up.set_data(x_up,y_up) \n",
    "    line_down.set_data(x_down,y_down)# update the data\n",
    "    time_step.set_text('Time step = {:}'.format(i)) \n",
    "    return line_up, line_down,\n",
    "\n",
    "ani = animation.FuncAnimation(fig, animate, init_func=init ,frames=n_steps,\n",
    "                              interval=1,blit=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save the animation as a gif"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ani.save('glide.gif',writer='imagemagick',fps=30)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}