Heat-Equation/main.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#1 Loading functions  and modules\n",
    "from fenics import *\n",
    "import matplotlib.pyplot as plt\n",
    "T = 0.1\n",
    "num_steps = 20\n",
    "dt = T/num_steps\n",
    "rho = 7500\n",
    "Cp = 500\n",
    "k = 50\n",
    "alpha = k/(rho*Cp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAACgCAYAAAAM9xtHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzsvXl0G9eZJX6LWEiCG0CCxMJ9AcEFZIE7JWs1rcWyLMmSZcuSvMiRJVmy40WxezoTZ+vjOHEnnenTm9uJx2c880t34nYcJUrS09NuO9pFiiLADaRIivtOgAAXECBB1u8P6j1XAQWJluUlCb5zdI4AfgCqLh7qvfru9+5lOI5DKEIRilCEIhRhX/QBhCIUoQhFKL4cEZoQQhGKUIQiFABCE0IoQhGKUITiRoQmhFCEIhShCAWA0IQQilCEIhShuBGhCSEUoQhFKEIBIDQhhCIUoQhFKG5EaEIIRShCEYpQAAhNCKEIRShCEYobIf2iD+CThFqt5jIyMr7owwhFKEIRij+aqK+vn+A4LnEluX9UE0JGRgauXLnyRR9GKEIRilD80QTDML0rzV1RyYhhmK0Mw7QzDNPJMMx/E/l7OMMwP7/x98sMw2Tw/vaXN55vZxhmy43nIhiGqWUYxsowTAvDMN9Z6QGHIhShCEUoPpu45YTAMIwEwD8AuBdAAYBHGIYp8Ev7CoBJjuNyAPwYwA9uvLYAwD4AhQC2AvjHG+/nBXA3x3EsADOArQzDVN+ZUwpFKEIRilDcTqzkDqESQCfHcdc5jpsH8K8Advrl7ATwv278/98A1DAMw9x4/l85jvNyHNcNoBNAJbccMzfyZTf+fWayqytVdP2y532Rn/1lz/siP/vLnvdFfvYXmffHcIxftlgJh5AMoJ/3eABAVbAcjuN8DMO4ACTceP6S32uTAXrnUQ8gB8A/cBx3+XZO4FbBcRy++93v0sc6nQ7Lc1VgDA0NfeI8jUYDiURyyzy1Wg25XH7LPKVSCYVCccu8qKgoxMXFieaNjY3B5/MBAGQyGRITxfkkjuMwPDxMH6/0nLVaLcLCxNcSX3ZsJicnMTc3B+Dm2Pi/553GJikpCVKp+M+Pn5eQkIDw8PBb5sXFxSEqKuqWeQqFAkqlUjRvZmYGU1NTAACpVIqkpCTRPP/3/FPBJjIyEiqVSjRvbm4Ok5OTAACJRAKNRiOa5/+edxKbb33rW0E/807FF0Yqcxy3CMDMMIwSwPsMw5g4jmv2z2MY5giAIwCQlpb2iT9nYWFB8JhcALVaLWJiYvjHI5qXkJCA+Pj4oO8/OjoKAIiJiYFWqw2aNzExAQCQy+VIS0sLOkicTiecTicAICsrK+gFdXZ2FrOzswCWyXaZTCaat7CwQAdUSkoKIiMjBX/jx0qxGRkZAQDEx8cjISFB9HOBj7GJjo6GVqsNes4EG5lMhvT09M8UG/J64ObYkAmVBMFGo9EgNjY2yBl/jI1KpYJarQ6aNzY2BmB58rrZRcNutwNYvkBnZGQEzXO5XHC5XABujo3b7Ybb7QYApKenCyZighmwfP4Em+TkZMFEvLi4KHhPgk1SUlLQiRj45NgoFAro9fpbYiORSJCZmfmpsZmbm6OLBX9svF4v/f/i4uJnho1SqbzpIuWzjpVMCIMAUnmPU248J5YzwDCMFEAcAPtKXstxnJNhmA+xzDEETAgcx70J4E0AKC8v/8T3WHK5HDKZDEVFRdDr9bBarejv78fo6ChiYmJgNpuRm5sLqVSKN954A9HR0TCbzbBarejq6oLdbkdUVBTMZjMKCgoQHh6Od999FyMjI9i6dSssFgva29sxPT0NhUIBs9mMoqIiREVF4YMPPsD58+fxxBNPwGKxoKWlBfPz87Db7WBZFizLQqlUwmKx4NSpU9i3bx96enrQ1NSE2dlZjI6Oori4GGazGUlJSejr68Pbb7+N7du3Y3p6GlarFU6nE8PDwygsLITZbEZKSgrcbjd++MMfYt26dYiIiIDFYsHY2BiGh4eRl5cHs9mMrKwshIWF4fvf/z7y8/ORmpoKq9WKvr4+jI6OUhyMRiOkUil+8pOfIDw8HGVlZbBYLOjq6oLD4aDYFBYWIjw8HO+99x4GBwexbds2WCwWtLW10ZUnwSY6Ohoffvghzpw5gyeffJJi4/V6MTExQbFRqVRobGzE+++/j4ceegj9/f1obGzE7OwsRkZGKDYajQaDg4P46U9/im3btsHtdsNiscDpdGJoaIhik5qaCo/Hg9dffx1r1qxBVFQULBYLRkdHMTw8DKPRCLPZjOzsbISFheH111+HwWBARkYGLBYL+vr6MDY2JsBGJpPhrbfegkQiQWVlJSwWCzo7OzE5OSnAJiIiAr/61a/Q3d2NHTt2UGxmZ2fhcrlgNptRXFyM6OhonDlzBh9++CEOHTqExsZGNDc3w+v1Ynx8nGITHx+P5uZmvPfee9i7dy8GBwfR2NiImZkZjIyMoKioCGazGVqtFsPDw3jzzTexdetWeL1eWCwWTE5OYnh4GAUFBTCbzUhLS8PCwgJee+01rF69GrGxsbBYLBgZGcHIyAiMRiNYlkVOTg7CwsLwox/9CBkZGcjOzobFYkFvb68Am7y8PMhkMrz99tvgOA6rVq2CxWJBR0cHJicn6W+lsLAQkZGROHXqFDo7O/HAAw/AYrHAZrPB7XbD5XKBZVkUFxcjJiYG586dwwcffIDHH38cLS0taG5uhsfjwdjYGFiWhdlsRnx8PGw2G37xi19gz549GBkZgdVqxczMDIaHhwXYjI2N4Y033sCWLVuwsLAAi8UCh8NBxw3LskhPT8fi4iJeffVVrFq1iv5mh4eH6bgh2EgkEvz4xz9GSkoKDAYDrFYrenp6KDYsyyI/Px8ymQzvvPMO5ufnsWbNGoqN0+lEVFQUWJaFyWRCZGQkTp8+jba2tk966butYG5Vx7pxgb8GoAbLF/M6APs5jmvh5ZwAUMRx3DGGYfYB2M1x3EMMwxQC+BmWeQg9gA8AGADEA1i4MRlEAvgPAD/gOO70zY6lvLycu5220+985ztYtWoVNm/eDGB5RWq1WtHY2IipqSlERESgqKgIdXV1MBqN2LdvHwBgamoKjY2NsFgssNvtkMlkyM/PR3NzM5RKJZ599lkAyyuL5uZmWCwWDA0NISwsDAaDASMjI3C5XPRWb2FhATabDRaLBd3d3QCWV7CLi4vo7+/Hc889B6VSicXFRXR2dsJqtaK9vR1LS0vQ6XRISEhAc3MzDh48iOzsbHAch97eXlgsFrS2tmJhYQHx8fEwGAy4fPkytm3bhoqKCnAch5GREVgsFjQ1NWFubg4xMTEoLi7G+fPnUVVVha1btwJYXnVZrVZYrVaKjclkwpUrV5CTk4MDBw4AAKanpyk2ExMTkEqlyM/PR2trK2JiYvDcc89RbFpaWmCxWDA4OAiGYWAwGDA2Ngan0ynApq2tDRaLBdevXwewvEoDgN7eXjzzzDNISEjA4uIiurq66ES8tLQErVaLxMRENDU1Yf/+/TAYDOA4Dn19fXSyIdjk5ubi0qVL2LJlC6qrl/sY+Ni43W5ER0ejuLgYFy5cQEVFBbZt2wYAcDgcFBuXy4Xw8HCYTCbU19cjKysLjz76KMWmqakJFosF4+PjkEqlyMvLQ3t7OyIiIvDiiy8CADweD8VmYGAADMMgJycHdrsdDocDr7zyCsLCwrCwsID29nY6ERNswsLC0N3djePHjyMxMRFLS0sCbBYXF6HRaKDVamG1WvHwww8jLy8PHMehv79fsEhRqVTIz8/HhQsXsGnTJqxevZpiQ34rbrcbUVFRKC4uxsWLF1FaWor7778fwHIpjmDjdDoRHh6OwsJCXL16FRkZGXj88ccBLJelCDZjY2OQSCTIy8tDR0cHZDIZvva1r1FsWltbYbFY0N/fT7FxOp0YHx/H17/+dchkMvh8PgE2HMchLS0NMpkMXV1dOHbsGDQaDZaWlnD9+nU6ERNsdDodLBYL9u7di4KCAnAch4GBAVgsFjQ3N2N+fh5KpRKFhYU4f/48ampqsGbNGgD
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#2 Create mesh and define function space\n",
    "nx = 0.008\n",
    "ny = 0.003\n",
    "mesh = RectangleMesh(Point(0,0),Point(nx,ny),30, 30,'left')\n",
    "V = FunctionSpace(mesh, 'Lagrange', 1) #Lagrange are triangular elements\n",
    "plot(mesh)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Boundary conditions\n",
    "u0 = Constant(100)\n",
    "def boundary(x, on_boundary):\n",
    "    return on_boundary\n",
    "\n",
    "bc = DirichletBC(V,u0, boundary)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "u_n = project(1, V)\n",
    "u = TrialFunction(V)\n",
    "v = TestFunction(V)\n",
    "f = Constant(0.0)\n",
    "F = u*v*dx + alpha*dt*dot(grad(u), grad(v))*dx-u_n*v*dx\n",
    "a, L = lhs(F), rhs(F)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "vtkfile = File('solution/solution.pvd')\n",
    "u = Function(V)\n",
    "t = 0\n",
    "for n in range(num_steps):\n",
    "    t += dt\n",
    "    #u0.t = t\n",
    "    solve(a == L, u, bc)\n",
    "    #c = plot(u,)\n",
    "    #plt.colorbar(c)\n",
    "    #plt.show()\n",
    "    ####\n",
    "    vtkfile << (u, t)\n",
    "    u_n.assign(u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1e-13"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1E-13"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "100.0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1e2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "100.0000000000001"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1E-13+1e2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Boundary conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "#1 Loading functions  and modules\n",
    "from fenics import *\n",
    "import matplotlib.pyplot as plt\n",
    "T = 4\n",
    "num_steps = 200\n",
    "dt = T/num_steps\n",
    "rho = 7500\n",
    "Cp = 500\n",
    "k = 50\n",
    "alpha = k/(rho*Cp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAACgCAYAAAAM9xtHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJztnXtwVOeZp58PXRAIgdAdIYQkJCQEUrcBc4mxjQ0GTEiYZDKBpTKTZL3lqUom2WS3KjXZmTJjT6UymXUlQ9Uk3nGFbCbeSYztxAlDsLmZYAyYq8+R0I07RoAQCIEkLlJL+vaP7j5uSS2pjyyd04L3qaKq+5zTRz9d6O/0873fe5TWGkEQBEEY43YAQRAEITqQAUEQBEEAZEAQBEEQAsiAIAiCIAAyIAiCIAgBZEAQBEEQABkQBEEQhAAyIAiCIAiADAiCIAhCgFi3A9ghLS1N5+XluR1DEARh1HD8+PEbWuv0SI4dVQNCXl4ex44dczuGIAjCqEEpdTHSYyNSRkqpVUqpOqXUGaXU34bZP1YptSWw/7BSKi9k3/cD2+uUUisD2xKUUkeUUqZSqkop9WKkgQVBEISRYdABQSkVA/wUeBYoBf6LUqq012HPAc1a60LgJ8CPAq8tBdYDs4FVwM8C52sHntZaewAvsEoptWh4viVBEARhKETyCWEBcEZrfU5r3QG8Dqztdcxa4N8Dj98ClimlVGD761rrdq31eeAMsED7aQscHxf4N2JtV6O1o6vksofksofkske05nKSSOYQpgKXQp7XAwv7O0Zr3amUug2kBrZ/2Ou1U8H65HEcKAR+qrU+PJRvYDC01rz00kvW8ylTpuAfq9znypUr1mPJNTijIVdWVhZjxkRH8d5oyJWZmUlMTIyLaT4hNFdGRgaxsdExxRrMtXHjxhH/Wq59x1rrLsCrlEoG3lZKzdFan+x9nFLqeeB5gNzcXNtfx+fz9Xh+9epVwP8fJCkpyX7wYaKzs7PH82CuzMxMJk6c6EYkALq7u3s8D+ZKT08nOTnZjUhA36u3YK7U1FRSUlLciBSWhoYGACZPnkxaWprLaT4hmGvSpElkZGS4nOYTrl27BkBSUhJZWVkup/mExsZGABITE6Pq4mOkiWRAuAxMC3meE9gW7ph6pVQsMAloiuS1WutbSqm9+OcY+gwIWutXgVcB5s+fb/szXXx8PHFxcZSVlZGdnY1pmly6dIlr166RlJSE1+tl5syZrlwNvPzyyxQUFJCfn49pmly8eJFr164xYcIEPB4PJSUlxMXFOZ5r06ZNTJkyhZKSEgzD4Pz581y/ft3KNWvWLOLj4x3P9corrzBx4kQ8Hg+GYXD27FmamppITEzE6/VSWlrK2LFjHc+1efNmYmJiWLBgAYZhcObMGZqbm61cs2fPJiEhwfFcr732Gvfv3+fxxx/HNE1OnTrF7du3rVxz5sxh3LhxjufasmULTU1NLF++HNM0qauro7W1lcTERDweD2VlZSQmJjqe6+233+bixYusWbMG0zSpqanhzp07tLS04PF4KC8vZ8KECY7n2rZtG7W1tY58rUjeBY8CRUqpfPxv5uuBDb2O2Qp8FTgEfAl4T2utlVJbgV8rpX4MZANFwBGlVDrgCwwG44BnCExEjxRjx45l3rx5zJs3j6amJgzDoKKigjfffJOEhATmzJmD1+slOzvb0auB+Ph4HnnkER555BFu3ryJaZqYpsnvfvc7xo4dy+zZs/F6veTk5DiaKy4ujvLycsrLy7l165aV6/e//z3bt2+ntLQUr9dLbm6uo7liY2OZM2cOc+bMoaWlhYqKCgzDYOvWrbzzzjvMmjULr9dLXl6eo7liYmIoLS2ltLSU1tZWKisrMQyDbdu28e6771JSUoLX6yU/P99RdRMTE0NJSQklJSXcuXPHyrV9+3Z27NhBcXExXq+XGTNmOJprzJgxzJw5k5kzZ3L37l1OnjyJaZrs2LGDXbt2MXPmTDweD0VFRY4qpTFjxlBYWEhhYSH379+3cu3atYvdu3dTVFSEx+Nx7SJypBn0OwrMCfwNsAOIAX6hta5SSr0EHNNabwU2A68ppc4AN/EPGgSOewOoBjqBb2qtu5RSU4B/D8wjjAHe0FpvG4lvMBypqaksW7aMp556ivPnz2OaJoZhcOzYMdLT062rAaeVUkpKCk899RRLly7lwoULmKZJZWUlJ06cIDU1FY/Hg8fjcVwpJScn8+STT/LEE0/w8ccfYxgG1dXVGIbB5MmTrVxOK6WJEyeyZMkSHnvsMS5fvoxhGJw8eZKKigomTZpk5XJaKSUlJfGZz3yGxYsXc/XqVQzDoLKykpMnT5KUlGTlclopJSYmsmjRIhYtWkRDQ4OVq7q6mgkTJlBeXo7H43FcKY0fP54FCxawYMECGhsbrYu12tpaxo8fT1lZGV6v13GllJCQwPz585k/fz43btywcp06dYpx48ZZF5EPklKKaIjTWm8Htvfa9kLI4/vAX/Tz2h8AP+i1rQJ4xG7Y4WbMmDHMmDGDGTNmcP/+faqqqjBNk927d7Nnzx4KCwtdUUpKKfLz88nPz+fZZ5+luroa0zR57733eO+995gxY4YrSkkpxfTp05k+fTrPPvsstbW1GIbBn/70J/70pz+Rn5/vilJSSpGTk0NOTg4rV66krq4OwzB4//33ef/998nNzXVFKSmlyM7OJjs7mxUrVnDq1CkMw+DAgQN88MEH5OTkuKaUsrKyWLVqFc888wynT5/GNE0+/PBDDh48SHZ2tmtKKSMjgxUrVrB8+XLOnDmDaZocO3aMw4cPk5WV5ZpSSktLY/ny5Tz99NOcO3cO0zQ5ceIER48eJSMjw1WlNJw8eJ95hkhCQkJUKqWxY8dGpVKKj4+PSqUUFxcXlUopNjZWlJINRCm5w+hL7ACilOwhSskeopTsIUrJOWRAGABRSvZziVKyl0uUkj1EKY0sMiBEiCgle4hSsocoJXuIUhoZoi/RKECUkj1EKdlDlJI9RotSCl4URbNSkgHhUyBKyX4uUUr2colSskc0K6XQi0jDMKJSKcmAMEyIUrKHKCV7iFKyRzQrpdCLyGhTSjIgjACRKiWniVQpOU2kSslpIlVKThOpUnKaSJWS00SqlJwmUqXkZBdWNZpavs6fP18P5Y5pL77ov//OokWLXOus2N7eTlVVFffu3euzz81cHR0dVFdXc+fOnT77Fi5c6NrEl8/no7a2lpaWlj77FixY4EqPJ/A3Jayrq+PWrVt99j366KOu9HgC6Orq4tSpU9y8ebPPvnnz5rnSSwn8uc6cOcONGzf67Js7d64rvZTA38Tx7NmzVhO7ULxeryu9lMDfxPH8+fNW88ZQ/v7v/35I7xNKqeNa6/kRHfugDwgdHR388Ic/7LPdzZa7XV1d/e4bM2aMaxNMksseozGXUsq1NtijMRe4917RO9cLL7wwpL8pOwPCA6+M4uPjiYmJISUlhfT0dOrq6ujq6iItLc1SN25dDfzoRz8iJiaG/Px8amtr6ezstNSNG1VKQX7yk5/Q3t5OSUkJ1dXV+Hy+HurGrfbcP/vZz7h58yZlZWVUV1fT0dHRQ9241Z578+bN1NfXM3fuXKqqqmhvb++hbtxqz/2rX/2K8+fP8+ijj1JZWcn9+/etjrVuVCkF2bJlC7W1tSxcuJDKykru3r1LYmKipW7cas/99ttvU1FRweLFi6moqODOnTuuVikF2bZtG6ZpOnKB8cAPCPBJbfCKFSu4d+8eJ0+exDAMdu7caU3keL1exyeYghNya9assaqUDMPoUaXk8XgoLi52VN2MGTOG4uJi/uzP/syqUjIMIyqqlIqKili7di3PPvssNTU1mKZpVSnl5eXh9Xpdac9dUFDA5z73OVatWkVtbS2maUZFldK0adNYvXp11FUpZWZm9qhSMgyDw4cPc+jQIbKzs61qIKeV0uTJk1mxYgXLli3j7NmzGIbB0aNHXa9Scurv5qEYEEIZN24cjz76KI8++iiNjY2YpklFRQV1dXWu1ywHq5RCJ5hOnz4tVUphiI+Pt650o61Kqay
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#2 Create mesh and define function space\n",
    "nx = 0.008\n",
    "ny = 0.003\n",
    "mesh = RectangleMesh(Point(0,0),Point(nx,ny),8, 8,'left')\n",
    "V = FunctionSpace(mesh, 'Lagrange', 1) #Lagrange are triangular elements\n",
    "plot(mesh)\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Boundary Conditions\n",
    "tol = 1E-14\n",
    "\n",
    "def BC1(x, on_boundary):\n",
    "    return on_boundary and abs(x[0]-nx) < tol\n",
    "\n",
    "def BC2(x, on_boundary):\n",
    "    return on_boundary and abs(x[0]-0) < tol\n",
    "\n",
    "bc1=DirichletBC(V,Constant(25),BC1)\n",
    "bc2=DirichletBC(V,Constant(800),BC2)\n",
    "bc=(bc1,bc2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "u_n = project(25, V)\n",
    "u = TrialFunction(V)\n",
    "v = TestFunction(V)\n",
    "f = Constant(0.0)\n",
    "F = u*v*dx + alpha*dt*dot(grad(u), grad(v))*dx-u_n*v*dx\n",
    "a, L = lhs(F), rhs(F)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "vtkfile = File('solution/solution2.pvd')\n",
    "u = Function(V)\n",
    "t = 0\n",
    "for n in range(num_steps):\n",
    "    t += dt\n",
    "    #u0.t = t\n",
    "    solve(a == L, u, bc)\n",
    "    vtkfile << (u, t)\n",
    "    #c = plot(u,)\n",
    "    #plt.colorbar(c)\n",
    "    #plt.show()\n",
    "    u_n.assign(u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}