{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Time-dependent Heat or Diffusion Equation\n", "\n", "#### Basic idea behind solving time-dependent (first order in time) problem in FEA\n", "\n", "A simplest version of time-dependent heat equation for temperature $T$ (or diffusion with concentration $c$) the following form:\n", "\n", "\\begin{align}\n", "{\\partial T(x,y,t) \\over \\partial t} &= k \\nabla^2 T(x,y,t) + s(x,y,t)\\;\\;\\;\\;\\mbox{ over the domain}\\\\\n", "T(x,y,t) &= T_D\\;\\;\\;\\;\\mbox{ on the boundary}\\\\ \n", "T(x,y,0) &= T_0(x, y)\\;\\;\\;\\;\\mbox{ initial condition}\n", "\\end{align}\n", "\n", "Here, we note that $\\nabla^2 = {\\partial \\over \\partial x^2} + {\\partial \\over \\partial y^2}$. The additional component that is present in the current problem is the initial condition, which defines what was the configuration of the system at $t = 0$, the initial time. Note that the boundary conditions could be of many variety as discussed in the earlier lectures. Here, however, we will only look at the simplest boundary condition. The question is how to solve the problem using FEA?\n", "\n", "We utilize a very simple trick in which the derivative in time is approximately discretized. One of easiest way to express the first derivative is as follows.\n", "\n", "$${\\partial \\over \\partial t}T(x,y,t) = k \\nabla^2 T(x,y,t) + s(x,y,t),$$ \n", "\n", "is approximately written as:\n", "$${T^{t + \\Delta t} - T^{t} \\over \\Delta t} = \\theta (k \\nabla^2 T^{t+\\Delta t} + s^{t+\\Delta t}) + (1-\\theta) (k \\nabla^2 T^{t} + s^{t}),$$\n", "where the super-script $t$ is the time and $\\theta$ is a parameter that can vary between $0$ and $1$. This equation can be re-written to give an equation of the form:\n", "$$T^{t + \\Delta t} - T^t - \\Delta t[\\theta (k \\nabla^2 T^{t+\\Delta t} + s^{t+\\Delta t}) - (1-\\theta) (k \\nabla^2 T^{t} + s^{t})] = 0.$$ \n", "What this means is that if we have somehow found $T^{t}$ then we get an $T^{t + \\Delta t}$ solving the above boundary value problem, subject to the boundary conditions. It is obvious that since the initial condition is known at $t = 0$, we can get started with the problem, and evaluate the value of $T$ at $\\Delta t$, and that way can proceed further. \n", "\n", "#### Special cases for different values of $\\theta$\n", "When $\\theta = 0$, this scheme is called as _explicit euler_ scheme. It means that by knowing the solution at $t = 0$ (initial condition), we can simply calculate what happens at $t + \\Delta t$, and so on. The explicit schemes are generally known to be unstable for larger time-steps (we will see that in the code below.) The convergence rate for this scheme is $\\Delta t^{1}$\n", "\n", "When $\\theta = 1$, the scheme is called as _implicit euler_. The scheme is known to be unconditionally stable, but the convergence rate is the same as that of explicit euler. \n", "\n", "When $\\theta = 0.5$, this scheme is called as [Crank-Nicholson](https://en.wikipedia.org/wiki/Crank–Nicolson_method) method. The convergence rate in many cases could now go as $\\Delta t^2$. Check this [link](http://www.mathematik.uni-dortmund.de/~kuzmin/cfdintro/lecture8.pdf) for additional details. \n", "\n", "To summarize, we use a __finite difference__ scheme for time-derivative and solve the resulting problem as a BVP using FEA.\n", "\n", "#### Method of manufactured solution for time-dependent problem\n", "\n", "We solve this problem using FEniCS in the following manner. Consider a unit square domain. If we choose the conductivity $k = 1$, the BVP becomes\n", "\n", "\n", "\\begin{align}\n", "{\\partial T(x,y,t) \\over \\partial t} &= \\nabla^2 T(x,y,t) + s(x,y,t)\\;\\;\\;\\;\\mbox{ over the domain}\\\\\n", "T(x,y,t) &= T_D\\;\\;\\;\\;\\mbox{ on the boundary}\\\\ \n", "T(x,y,0) &= T_0(x, y)\\;\\;\\;\\;\\mbox{ initial condition}\n", "\\end{align}\n", "\n", "and with a special case:\n", "\n", "\\begin{align}\n", "T(x,y,t) &= 1 + x^2 + \\alpha y^2 + \\beta t\\\\\n", "s(x,y,t) &= \\beta - 2 - 2 \\alpha\n", "\\end{align}\n", "\n", "which clear satisfies the strong form of the differential equation, with the initial condition corresponding to $t = 0$ and the boundary conditions as per $t$ and $x,y$ on the boundaries. We now, solve this problem using FEniCS as follows:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### FEniCS solution\n", "The problem is taken from the FEniCS tutorial. The [program](https://github.com/hplgit/fenics-tutorial/blob/master/pub/python/vol1/ft03_heat.py) is kept on github. \n", "\n", "Below, we demonstrate the steps involved. First, we load all the libraries and define the relevant parameters required to solve. Please note that we use the variable $u$ for trial function and $v$ for test function. " ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "\"\"\"\n", "FEniCS tutorial demo program: Heat equation with Dirichlet conditions.\n", "Test problem is chosen to give an exact solution at all nodes of the mesh.\n", " u'= Laplace(u) + f in the unit square\n", " u = u_D on the boundary\n", " u = u_0 at t = 0\n", " u = 1 + x^2 + alpha*y^2 + \\beta*t\n", " f = beta - 2 - 2*alpha\n", "\"\"\"\n", "\n", "from __future__ import print_function\n", "from dolfin import *\n", "import numpy as np\n", "%matplotlib inline\n", "\n", "T = 2.0 # final time\n", "num_steps = 10 # number of time steps\n", "dt = T / num_steps # time step size\n", "alpha = 3 # parameter alpha\n", "beta = 1.2 # parameter beta" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, we define geometry and create the mesh. This is followed with creating a function space for the test (weight) and trial (temperature) functions. " ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "# Create mesh and define function space\n", "nx = ny = 8\n", "mesh = UnitSquareMesh(nx, ny)\n", "V = FunctionSpace(mesh, 'P', 1)\n", "\n", "# Define boundary condition. Note that the parameters can be passed as arguments to Expression\n", "u_D = Expression('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t',\n", " degree=2, alpha=alpha, beta=beta, t=0)\n", "\n", "def boundary(x, on_boundary):\n", " return on_boundary\n", "\n", "bc = DirichletBC(V, u_D, boundary)\n", "\n", "# Define initial value\n", "u_n = interpolate(u_D, V)\n", "#u_n = project(u_D, V)\n", "\n", "# Define variational problem\n", "u = TrialFunction(V)\n", "v = TestFunction(V)\n", "f = Constant(beta - 2 - 2*alpha)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The weak form of the current problem will be:\n", "$$T^{t + \\Delta t} - T^t - \\Delta t[\\theta (k \\nabla^2 T^{t+\\Delta t} + s^{t+\\Delta t}) - (1-\\theta) (k \\nabla^2 T^{t} + s^{t})] = 0.$$ \n", "\n", "Note that we define the initial value for temperature in variable `u_n`. This variable corresponds to $T^{t}$ in the theory above and the variable `u` corresponds to the solution at the current time $t + \\Delta t$. Also, `f` corresponds to the source term $s$.\n", "\n", "$$\\int_\\Omega T^{t+\\Delta t} w d\\Omega - \\int_\\Omega T^{t} w d\\Omega + \\theta \\left(\\int_\\Omega k\\nabla T^{t+\\Delta t} \\cdot \\nabla w d\\Omega - \\int_{\\Omega} s^{t + \\Delta t}w d\\Omega\\right)\\Delta t + (1-\\theta) \\left(\\int_\\Omega k\\nabla T^{t} \\cdot \\nabla w d\\Omega - \\int_{\\Omega} s^{t}w d\\Omega\\right)\\Delta t = 0\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, define the weak form in FEniCS. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "theta = 0.5\n", "F = u*v*dx + theta*dt*inner(grad(u), grad(v))*dx + (1. - theta)*dt*inner(grad(u_n), grad(v))*dx - (u_n + dt*f)*v*dx\n", "#F = u*v*dx + dt*inner(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx\n", "#F = u*v*dx + dt*inner(grad(u_n), grad(v))*dx - (u_n + dt*f)*v*dx\n", "a, L = lhs(F), rhs(F)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that `F` corresponds to the expression for the weak form, for particular value of $\\theta$, which equals $0.5$ correponding to Crank-Nicolson scheme. Note that in FEniCS `lhs(F)` correspond to all bi-linear term (both `u, v` or `grad(u), grad(v)`. In the current problem:\n", "\n", "`lhs(F) = a = dt*inner(grad(u), grad(v))*dx + u*v*dx`\n", "\n", "everything else will be stored in `rhs(F) = L`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, we update the solution in time. The term `u_D.t = t` ensures that the boundary conditions are updated corresponding to time `t`. The command `solve(a == L, u, bc)` will save the solution in `u`. The solution that is obtained is now used to update the previous solution `u_n`. The solution then continues as we discussed above. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Time-stepping\n", "file = File(\"./Diffusion/output.pvd\", \"compressed\")\n", "\n", "u = Function(V)\n", "t = 0\n", "for n in range(num_steps):\n", "\n", " # Update current time\n", " t += dt\n", " u_D.t = t\n", "\n", " # Compute solution\n", " solve(a == L, u, bc)\n", " \n", "\n", " # Plot solution\n", " plot(u)\n", " #plt.pause(0.5)\n", "\n", " # Compute error at vertices\n", " u_e = interpolate(u_D, V)\n", " error = np.abs(u_e.compute_vertex_values(mesh) - u.compute_vertex_values(mesh)).max()\n", " print('t = %.2f: error = %.3g' % (t, error))\n", "\n", " # Update previous solution\n", " u_n.assign(u)\n", " file << (u, t)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The output is saved in [Paraview](https://www.paraview.org). You can download it. The solution from each step is stored in `.vtu` file. The full code is given below. More details of the solution could be obtained from the [FEniCS tutorial](https://www.springer.com/gp/book/9783319524610) that is freely available. Look at Chapter 3 of that book.\n", "\n", "The full code is given below. " ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Calling FFC just-in-time (JIT) compiler, this may take some time.\n", "Calling FFC just-in-time (JIT) compiler, this may take some time.\n", "t = 0.15: error = 3.55e-15\n", "t = 0.30: error = 2.66e-15\n", "t = 0.45: error = 4.44e-15\n", "t = 0.60: error = 3.55e-15\n", "t = 0.75: error = 4e-15\n", "t = 0.90: error = 3.55e-15\n", "t = 1.05: error = 4e-15\n", "t = 1.20: error = 4.44e-15\n", "t = 1.35: error = 4.44e-15\n", "t = 1.50: error = 4.44e-15\n", "t = 1.65: error = 6.22e-15\n", "t = 1.80: error = 5.33e-15\n", "t = 1.95: error = 4.44e-15\n", "t = 2.10: error = 5.33e-15\n", "t = 2.25: error = 6.22e-15\n", "t = 2.40: error = 7.99e-15\n", "t = 2.55: error = 7.11e-15\n", "t = 2.70: error = 7.11e-15\n", "t = 2.85: error = 7.11e-15\n", "t = 3.00: error = 6.22e-15\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQYAAAD8CAYAAACVSwr3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAer0lEQVR4nO2dbaylVXXHf2vunRcICAho7DBGavCFNpIogjF9QYwVaJqJqUkRU1OimZCI6UdIP2gTv2jsBzWikwmlxLRx0kRjp80oMW0UE6QOTRF5CWSKLVwxIYCFijNzZ+aufrj3js+ce16eZ++91t77nP1LTjLn3ueeZ9079/zuf629n3NEVWk0Go0u23IX0Gg0yqOJodFobKGJodFobKGJodFobKGJodFobKGJodFobGGmGETkHhF5XkQenfB5EZEvi8hREXlERN6ZvsxGo+FJn8RwL3DDlM/fCFyxcdsHfC2+rEajkZOZYlDV+4GXphyyF/i6rvMgcKGIvCFVgY1Gw5/lBI+xG3i2c39l42O/GD1QRPaxnipYPmf5XRe+6YIEp6+fJVnzPyc250z5vcQ+1jKnXc+9zPBdxEsy8HiGfcF/PrL6gqpeOuwsacQwrtKxPyFVPQAcAHjt2y/V6//2TxOcvh8XbD/mdq5UvGb7cb9zLaX/+bx2+VfJHuuCpV9Hff3FgbVcuC3svBdvG/Z/d8HAZYALtu3sddy5v/Xf/zPskddJIYYVYE/n/mXAcwkeNykvnzwn+GtzSeWVk7tmHpNKHq+cnvzzCZXGS6fO2/KxUFm8fPrcs+4PFcWLI7X0FcX/rp193r6ieHHtN/93fSTxcieg9JHEy2snOsf3k8QQUojhEHC7iBwErgVeVtUtbUTNhEjFSyaT5JEybYyTRm5ZlCCKeZbETDGIyDeA64BLRGQF+AywHUBV9wOHgZuAo8CvgVuTVFY5fWRiKY9paSOFNEqTRUpRWKaJWiQxUwyq+pEZn1fgk0NPfFq38auTO4Z+WTTnbV91P+ckcsnDKmVYysJTFF5pwlMSQ0nRSlRFjIxySMVTHhYpI5UsYlNFKlHUKomhLJwYYgiRiodMJskjZdoYJ43SZOEhipA0MbTlKEESTQzG9JGJlTympY0U0ihNFrlFYZEmckkimxjWVPi/k+mXWfpy/vbw/is1OeRhlTIsZZFLFPMgiaEsbGKIkVIOqcySRypxWAgjlSxyiSIkTcS2HKklMZSFFUMMQ6XiIZJJ4ihVGClkUYIorNJEbknkayXWhGOr2zlnx8lcJbjRRyRW8piWNFJIY5wwUskipyjmQRIxZE8Mx1a3Zz1/KWLKIQ+rlJFKFrGpIkYUIWkituVILYkYsoshNzFi8pbKNHmklIaFMKxkkUMUVmmiJEksvBhiGCIVa4l4SGOcMFLKohZRxEoCZositySyLlcePzn5ibVrexkRPxV9JGIlj0nSSCGMlLJIkSq8ReExwIyRRCjFJoZp0rCgBBF5y8MqZVjKIpcovNJEakmEUqwYvIkRkadUpsnDQxqhwkgli1yiiJUEDB9g5pREE0MC+krFWiAe0hgnjJyyyCEKr9lETklkE4OqcGJ1+Ol37jhlUI0PfQRiJY9J0kghjJJkkVIUXmnCUhKhVJcYQmQyhNzi8ZaHVcqwlIWXKLzShKUkQqlODNbEiMdLKl6rOalTRipZpBQF9JeFR5ooRRJNDAnpIxVreXhIY5wwcsoiRhRwtiyaJNbJOGOAU6tLyR5veUf4ewh4klMek6SRQhiWsvAURUjrEdtyWEsihLlJDCklA3lF4y0Pq5SRShaliMIqTVhLIoS5EUNqYkTjIZVp8rCWRm5Z5BJFrCRgtihKkUQTgwF9pGIpD2tplCaLHKIIXe0YkiZySiKfGFRYO5k2/nfZtr3smcMseViJY5I0YoVhKYucicIrTVhKIoS5TQwx0ilBKt7iGCeM1LIIFcVoovAUhbckLGYSIcytGGIIkYq3TCaJI6UwUssiVaqISRNwtihKk0To4DK1IJoYEtFHJh7yGCeMmmQxVBQp00SpksiRIvKJYQ04MeaF73dGvOZ14UyTh6U0rNPFqCxqFUUOSVimiBjKSwzjZGFFQRKalTgsxGElDCtReLYdXpKwbjVCKU8MnoRIKJNMPNNG6nYkVfuRK02ULgkLQSy2GEIYIhMniUySRkphWMsiVhReacJbErlSRBODJbMkYiyOccKwlEVOUaRME16SKDlFZN3gJKs+8wTdUc4s4SwyiMMyXZQqilIl4ZUiQliIxBAioCJkMk0ciaVhIYyU7UdXFJ5pwlsSFikihIUQQwjFy2SSNByEkVIWIaLI1XZ4SMJ6YNmXJoaEzJKJizgchDEqi5JE4TXELFkSKQTRSwwicgPwJWAJuFtVPzfy+QuAvwfeuPGYf6Oqfzf1QRW2rcrUQ9Z2aJ/yqiGrOAw3k5UkihyzCU9JhFz5GcJMMYjIEnAX8AFgBTgiIodU9fHOYZ8EHlfVPxGRS4EnReQfVDXq3VFniSMFJcnHXRxG6SJl+5FSFKVLwiJFhNInMVwDHFXVpwFE5CCwF+iKQYHzRUSA84CXgCpe5z1EPrlkMk0cSaVhIIxUqaIrCs804SEJyxQxlD5i2A0827m/Alw7csxXgEPAc8D5wJ+p6pbfIhHZB+wDWL7oopB6i2CITLwkMk4a5gkjsyhypQlPSeRKEX3EMO5ZMPrb/kHgYeB64M3A90Tkh6r6yllfpHoAOACwc8+ecjK8IdMkYi2NSQkjmTAKFkUNkrBMERAniT5iWAH2dO5fxnoy6HIr8DlVVeCoiPwMeBvw40kPKj2GjykpaZawSa7hq1m6SDjgjBVFTJrwkoRXqxFCHzEcAa4QkcuBnwM3A7eMHPMM8H7ghyLyeuCtwNPBVRlQ0yxhE09xmKWLRKmiKwrPNFG6JGLelXwaM8WgqqdE5HbgPtaXK+9R1cdE5LaNz+8HPgvcKyI/Zb31uENVXzCp2JHSZTKpPmthRMkigShypQlvSVikiL702segqoeBwyMf29/593PAHyWrqmJK2JsxrgZLWSQVBQyWRUpReEmi9BSR8SIq2DZgNrIW3i4VRS5xWKYL01QRmShKlYRXigilmi3RQyQylJKk472KYZUukqWKyNYjJk14S8IiRYRSjRgsCZFODpl4ScMiXZiIwjFNeEjCemA5hCaGQPrKxEsgHoPI0XOkFAUEyCJhmpgXSaQSRBODMbMEYi0Oy0Fk6seOThURaaIGSXimiGxikDVYMpwbjHK6oDlClxzisEwXVqnCM03kkIR1qzGUhUkMQyRUkkQmicNDGKlFkU0SEJwmYiVhIQgYliJCWBgxDKEGiYwTRmpZpG4VUokntyQsU4RlmzGEJoZI+kjESx7esohNFClEkUMSJaaI1IJoYnBgkjw8hGEpi9StR6x0kknCKUWUMIeYRDU7H/tS0malWUxLG5bSGP25lyiKrJJwShGhbUbI7sqhzF1iCJVNaULxTBlWqSJV21GEJIwHlkNSxFBBhDB3Ygillus2vISRWhap0kQ2STgNLEMEAekl0cQQQGm7HmG8MCxlEfu9pRBFdkkYpgjr5c5ZNDEY4rkHYRyjskgpitRzitgneSpJlJYicgki385HHb7zsaSNRzFMSxyW0rBMFSlFkVMSOVLEEEHA8DYjhKoSQ+ot1CWKxlsaVrJI1XaUIAnrFGE5hwilKjGkpoYdjl28WpPULUiqNJFLEl4poiRBLLQYhlCyRKz2JWxSoiiqksTAFGE9qOxDE4MBuTYubWK9NdpSFDVKopQUkVIQxex8LG2DkRW5tkdbyiKlKGLTRCpJuKUIQ0HEUExisHxNxy6lCihHyrBqQbrfS840ESMJtxQR2GbEvKN4H4oRgxe17HDsUutuR7BLEzkkUVKKsBbEwolhCENTjLdIvHc7QjmiyCEJ9xRhMKjsSxNDQnJtXOpiLQtLUdQoCdMUYbzcOY0mBidybo+uZcdjijThLQmXFJFBEHlfDPaE3/s8nt7p987aQyhle7SFKHKniRSS8EoRgwQBg9uMoSxMYhgioVIk4ikNi1RhlSY8JeGVIiznECEsjBiGUINEPF7fscQdj5BfEpYDy1IE0cQQySyJeIpj0XY8QjpJWKeI7tcNEQTYzCFm0cRgzDRxeEjD8jqKknY8jtYztJYcKcKkzUgkiKxbolNdRl3ClY8hTJKGpTC8tkbnThMpJOE1izAVRCBzkRhiBVOaWLxTxjzveBytZUgdXinCfA4RwFyIIZaS3jRmFuOkYS2LFIkiVZoI+Ws+rg7PFDG0zTBbyRhAr7whIjeIyJMiclRE7pxwzHUi8rCIPCYiP0hbZn6WVmffstV2Qs/cLNi2evYtlhQ/s9h6QmsIOee2VRn7dn+pjpfVbVveKTyWmYlBRJaAu4APACvAERE5pKqPd465EPgqcIOqPiMir0taZSXkfMepMzWMyKH0NJFyS7RnivBIEKYziBn0aSWuAY6q6tMAInIQ2As83jnmFuBbqvoMgKo+H13ZHOFxsdPEcxuLoiRJpJpFzIMgYuXQRwy7gWc791eAa0eOeQuwXUS+D5wPfElVvz76QCKyD9gHsOPciwZH31J2JKYgV7qwFIWFJHKkiHkQRGxr0UcM435zRitbBt4FvB84B/iRiDyoqk+d9UWqB4ADAOe9ds/ghjhVD12yYLzThZUoUkkiZ4ooWRCxbyA8iz5iWAH2dO5fBjw35pgXVPVV4FURuR+4CniKAukjmJLk4ZkuLESRWhLeKaJEQYSseAyhT944AlwhIpeLyA7gZuDQyDH/BPy+iCyLyLmstxpPpC3Vl+6kf9ItNx6rIqm/3xSrG7Hfc2gNMasYQ85luYLRl5mJQVVPicjtwH3AEnCPqj4mIrdtfH6/qj4hIt8FHgHWgLtV9dHk1RZGSddJgO1b0sHZ328pSSJHivDYC2E9oJxFrw1OqnoYODzysf0j978AfCFJVXOC12akiedPuEV5y2MXJokcF1PFCKLveXIJIuMLteSP4jnIJQvLNGEliRwpwksQQ88Ts806hKxbopdP2O317nJqZ9pdYanJcTGVlShSSqK2C6lKFEQoC3GtRB8BlSgPz3RhIQoLSdRwjUSIkEIEYSmHhRBDH2qRh5csUosilSRy7G7sntd6qdNyQDmEJoYBlCoPj+sjUg4yU0six6XYHoLI2V40MSRmkjw8hWEtitIkUcs1EkPPlXP+0MTgxDhheMnCUhSlSqIGQVjOH9aPDxdE3ne7jliVWCug348llyxSLy+eedyCJJFLEEPOaTl/WD8+/P+22sQQIxUoVyzerYhVmrCQRIwgwOcaie45SxFECNWKIZY+YilJHl7pwkIUqSSRK0WULAgrOSysGPpQujxGZVGDKFJLwjNFlCgIq/TQxBDJNHl4S8NbFCVIIkYQUO5FVJvnydVeNDEYMkkaXsKwFkVJksjRZnjshcg1f8h6EdXS8d+8TffpXeHvzFsb44ThIQtLUZQoidLnEENmBKGCCKWYxNCVhAWliyeHLKxEYSEJ7xRR4hbokMcPpRgxWNNHPKXJw7sV6YqiNEnEPiFCU8SiCmJhxNCHWuQxKgwLUVikiRSSyDWLKFkQFnJoYhjINHnkkkZXFLWkidjVBMiTIjy2QG+eJ2d6aGJIyDhpeMvCUxLzJAiPFqOm9NDEYMyoLDxFYS0JC0FAfJvhuZIRsvuwhvTQxOBMLlFYSqK0NsNbEF7pYejjx8gh69WVSyfClyhP78w/BExBjvbDQxLzIAiv9sJ65SKEahPDPEvFM1VYSaI0QZTcXpR08dQm1YohhmlSKVEaXqKwkETKNiNGEIueHoaykGKYRg3S8BCFpSRSCMK7vfBMD1aDySE0MQygVGlYX3OSWhIpBJFr/uAxnCyhtWhiSMQ4aeSQhZckUgoCwiWRUxDW7UW77HpOGZWFtygsJVFaisgxoPRoL3KlhyYGR3KKwkMStQvCOz2UPJjM+noMcvxUr2N113z6K5corCSxyILwmj2EvnL1UKp4xvUVyCxKF0yOOYWFJFK2GbkF4TWcLOVy603KfqYkZppgSpWGZ6qwlEStgvBqL0q53HqTMp8NGahFGl6iSC2JRRREzemhnN/4gilZGl1RWEsipSAgThI5BWHdXnhfMDWOJoZIJkkjhzCsJVFiikghiBLbi9ytRa+fpojcICJPishREblzynHvFpHTIvLhdCXWiRw/ddbNm6UTp8/cTB7/+OlkL+C77cRa9FsOLp9Ym/j2frNYOqFb3lin19cFTP2HnmvIysLSatwVlV1m/lkTkSXgLuADwApwREQOqerjY477PHBfmtLmi1E5eCYKyyRRWpsRkyBC2ot5TQ99fnrXAEdV9WlVXQUOAnvHHPcp4JvA873OrMq246tn3RaJXInCKklsJohSUkRoegCC08PQv9ZDz+OZHvr82doNPNu5vwJc2z1ARHYDHwKuB9496YFEZB+wD2DX8mu2fN5CDmu7jBd8E5EjUVglCYsUEZIgvNMDDP9rbT2YDKXPb9+4ikdV90XgDlU9LTL5G1TVA8ABgAvOecNwLQcwSzalimNcirCUhYUkUgsiV3vhtbRp1VqE0Oc3bQXY07l/GfDcyDFXAwc3pHAJcJOInFLVbyep0pBp4ihNGl6pIrUkUq1mxK5ihApiEdNDn9+sI8AVInI58HPgZuCW7gGqevnmv0XkXuBfapDCLEqXRlcUtUliUQRRa3qY+dukqqdE5HbWVxuWgHtU9TERuW3j8/vTllQHk6SRSxiekmiCsG8vQrZUp5RDr98gVT0MHB752FghqOpfxJdVL+OE4S0La0lYCALCJZFDEB7tRc7Wwv691xtZl2Qtl0NTL3vGLnfmWOL02BhltSFq6uOkeZhGX3JJojZBxNDkMLCYcY8R/xCBrK3B8RPrtwUlhyQsN1WlFESNchgqCGs51P9CLTFy2LUzXR0Z6crBayaxKYfUc4ilE6eLGFDGbo7y2PcQsqQZ8wbAfSlDDDEMkUolEvGWhIUgUg8oY1cvmhyGUb8YhjBOIoXLwlMSJQsiRXrw3FZduxza8HFzztG9FYrXTMJiBlHC7CFm1cJjKGl9EdYQmhjGUYEsahREquFkrsHkIslhsVqJGAptQzxajdQtRor2ImdrsQhtRUsMMRSWLKxThEWCiH6MDK2Fx3Jm7uTQxJCaAiRRkyBStBc5WouYl5HrfWxGOTQxWJJZEh6CSEWsIFJspx5K6OtM1iCHJgYvMkrCUhCltRfecgD7oWQOOTQx5CCTJGoRRIr0EEqTwzr5ViXWFI4dn37MObt8aslJVw5OqxybcrBYxZDjp5KuXoSuXMSsWoRupbZesQh9RagQyl6unCWOadQoFWdJWC11plzejL3uImY7dciSZmlyCKVsMcQwTSo1SCOTJEoUROy+h9j0sIhyWMwZw7HjW28l4ziTsJhDpJw9RH194OzBc69D72ONr9Sf38QwlHFyKDFZOCWJ1AkiZXrI0VqEJgcY/pJxJSSHJoZplC6LTUksmCBytRZe26hLkMNithIxlNiGOLQZqVuMFO1FjtbCaxt17raiiSEFpYjCYRaRUhAp9j7k2vMw73JoYrCgJEkYkTo9pBBE8Nc2OWyhicGa3JIwFERp7UWOay3mVQ5NDJ7klEQlgohNDzlai3mUQxNDLnJJwlgQqcidHoZSohxiBJHxWok1NPBJISUtGaag+3Pw+t6MljpTLm/GLm3GXmvhsZxpuZQZQ5X7GIYKpSqReEvCUBAp9z40OfhSVjVG6LHjZ25V4VmvQXuR4706x1FDWzEEy1eH3mQhxNClOkl4ziEKlkPOmUMIQ+VgOW8IYeHE0KUqSXgJwmA4WbscYl9Tsi+hryNpwUKLoUs1kvAUREIWUQ41txRNDGOoQhJNDlVg3VJY0cQwg6Il4ZEeErcWJcih9JZiCFapoZcYROQGEXlSRI6KyJ1jPv9REXlk4/aAiFyVvtT8FCuJytJDCasVJbcUJQwiZ4pBRJaAu4AbgSuBj4jIlSOH/Qz4Q1V9B/BZ4EDqQkujOEksoBxqWqmoraXokxiuAY6q6tOqugocBPZ2D1DVB1T1lxt3HwQuS1tm2RQjCa/WIhG1ysHzasy+5HiLut3As537Kxsfm8THge+M+4SI7BORh0TkoVUt5C9tYooRhCVzJocQPOYNOVNDHzGM25g9Vk8i8j7WxXDHuM+r6gFVvVpVr94hFW1TDiC7IJocelPyvGEI3u9duQLs6dy/DHhu9CAReQdwN7BXVV9MU179NDn0I/dA0ntnZF9yDSL7iOEIcIWIXC4iO4CbgUPdA0TkjcC3gD9X1aeSVDZHZE0P1nOHguRQy7yhhkHkTDGo6ingduA+4AngH1X1MRG5TURu2zjs08DFwFdF5GERecis4orJLggr5kgOIcxjS9HrWlZVPQwcHvnY/s6/PwF8IrqaBWFTDu6Xgx87bncp9/ETyS7djr1kO9dl2pZ4X5rddj5mJEuCsGwt5iQ5lNpSDCE2NWQTg6qyduz4ltsikqW9aHJIjvUSpucgsrhXcOojh201vSJTT7K0F1atRUFtRSheLcXQV3zyaimKE0Mf5lke7oKoQA4xeM8bvF4OzpqyqklI7W2K6/zBau6QqK1oLcXZeCxfzq0YxlGjJNwFkZrK5eC18cn6RV2GslBi6FKbJJoc6pJD7alhYcXQpRZJuKUHi9aiEDl4MlQOJe2IbGIYoQZJuAoiJYZvstsX75bCWg5WNDFMoXRJLKocamoprLFKDU0MPSlVEi7pocnhDB7zhhJSQxNDAKUKwpQ5lMO8YJEamhgiKE0O5szZ9+u5v6HEV5ieRhNDJAsnh4YLuduJJoYElCKH7K81OZQCVilCKHEImZomhkSUIgdzCvs+a9subUXqOUMTQ0JKkEN1qWGBqGnO0MSQmBLkUBWVthPzThPDHFLd0mUkOZYtPeYMOQeQTQwGtNTQqJ0mhkZ+MrcTbQC5lSYGI3KnhkVrJ2qhlgFkE4MhueWwSMzrnCEXTQxzTFVLl211Yiy5BpBNDMbMdWqYo+9tXuYMqWhimHOqSg2NaFINIJsYHJjr1JCSyHailsuwaxhANjE4MbdymNfvqyfzOoBsYlgAWjvRj1LnDDkGkE0MjsxtakjJgrQTpdPEsCCYpoYmvKJIMYBsYnCmpYb5o8Q3pImliSEDueRQzayhbXbKThNDIw0FSSdmztAGkBvncz1b4wytpWiUTC8xiMgNIvKkiBwVkTvHfF5E5Msbn39ERN6ZvtRGCqppJyqjtP0MsQPImWIQkSXgLuBG4ErgIyJy5chhNwJXbNz2AV+LqmpBmLvUkPL7WYA5Q8kDyD6J4RrgqKo+raqrwEFg78gxe4Gv6zoPAheKyBsS1zqX5JDDIqSGeZwzeLLc45jdwLOd+yvAtT2O2Q38onuQiOxjPVEAnPjeyW88OqjavFwCvGDyyCeTP+LsWl9Jfs4Y7H626ampVoC3hnxRHzHImI9pwDGo6gHgAICIPKSqV/c4fxHUVG9NtUJd9dZUK6zXG/J1fVqJFWBP5/5lwHMBxzQajUroI4YjwBUicrmI7ABuBg6NHHMI+NjG6sR7gJdV9RejD9RoNOpgZiuhqqdE5HbgPmAJuEdVHxOR2zY+vx84DNwEHAV+Ddza49wHgqvOQ0311lQr1FVvTbVCYL2iumUU0Gg0Fpy287HRaGyhiaHRaGzBXAw1bafuUetHN2p8REQeEJGrctTZqWdqvZ3j3i0ip0Xkw571jdQws1YRuU5EHhaRx0TkB941jtQy63fhAhH5ZxH5yUa9feZqJojIPSLyvIiM3RcU9BxTVbMb68PK/wJ+G9gB/AS4cuSYm4DvsL4X4j3Av1vWFFnre4GLNv59Y65a+9bbOe7fWB8Qf7jUWoELgceBN27cf13JP1vgr4DPb/z7UuAlYEemev8AeCfw6ITPD36OWSeGmrZTz6xVVR9Q1V9u3H2Q9f0auejzswX4FPBN4HnP4kboU+stwLdU9RkAVS29XgXOFxEBzmNdDFn2Uqvq/Rvnn8Tg55i1GCZtlR56jAdD6/g46xbOxcx6RWQ38CFgv2Nd4+jzs30LcJGIfF9E/kNEPuZW3Vb61PsV4O2sb+T7KfCXqprvfeunM/g51mdLdAzJtlM70LsOEXkf62L4PdOKptOn3i8Cd6jq6fU/bNnoU+sy8C7g/cA5wI9E5EFVfcq6uDH0qfeDwMPA9cCbge+JyA9VtayrUNYZ/ByzFkNN26l71SEi7wDuBm5U1RedahtHn3qvBg5uSOES4CYROaWq3/Yp8Qx9fw9eUNVXgVdF5H7gKiCHGPrUeyvwOV1v4o+KyM+AtwE/9ilxEMOfY8ZDkWXgaeByfjPE+Z2RY/6YswcjP840wOlT6xtZ39353hw1Dq135Ph7yTd87POzfTvwrxvHngs8CvxuwfV+DfjrjX+/Hvg5cEnG34c3MXn4OPg5ZpoY1G47da5aPw1cDHx146/wKc10pV3PeougT62q+oSIfBd4BFgD7lbVLJfl9/zZfha4V0R+yvoT7g5VzXI5toh8A7gOuEREVoDPANs7tQ5+jrUt0Y1GYwtt52Oj0dhCE0Oj0dhCE0Oj0dhCE0Oj0dhCE0Oj0dhCE0Oj0dhCE0Oj0djC/wNTcJjmuRIs8wAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "\"\"\"\n", "FEniCS tutorial demo program: Heat equation with Dirichlet conditions.\n", "Test problem is chosen to give an exact solution at all nodes of the mesh.\n", " u'= Laplace(u) + f in the unit square\n", " u = u_D on the boundary\n", " u = u_0 at t = 0\n", " u = 1 + x^2 + alpha*y^2 + \\beta*t\n", " f = beta - 2 - 2*alpha\n", "\"\"\"\n", "\n", "from __future__ import print_function\n", "from dolfin import *\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "\n", "T = 3.0 # final time\n", "num_steps = 20 # number of time steps\n", "dt = T / num_steps # time step size\n", "alpha = 3 # parameter alpha\n", "beta = 1.2 # parameter beta\n", "\n", "# Create mesh and define function space\n", "nx = ny = 8\n", "mesh = UnitSquareMesh(nx, ny)\n", "V = FunctionSpace(mesh, 'P', 1)\n", "\n", "# Define boundary condition\n", "u_D = Expression('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t',\n", " degree=2, alpha=alpha, beta=beta, t=0)\n", "\n", "def boundary(x, on_boundary):\n", " return on_boundary\n", "\n", "bc = DirichletBC(V, u_D, boundary)\n", "\n", "# Define initial value\n", "u_n = interpolate(u_D, V)\n", "#u_n = project(u_D, V)\n", "\n", "# Define variational problem\n", "u = TrialFunction(V)\n", "v = TestFunction(V)\n", "f = Constant(beta - 2 - 2*alpha)\n", "\n", "theta = 0.5\n", "F = u*v*dx + theta*dt*inner(grad(u), grad(v))*dx + (1. - theta)*dt*inner(grad(u_n), grad(v))*dx - (u_n + dt*f)*v*dx\n", "#F = u*v*dx + dt*inner(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx\n", "#F = u*v*dx + dt*inner(grad(u_n), grad(v))*dx - (u_n + dt*f)*v*dx\n", "\n", "\n", "a, L = lhs(F), rhs(F)\n", "#a = dt*inner(grad(u), grad(v))*dx + u*v*dx\n", "#L = (u_n + dt*f)*v*dx\n", "\n", "#a = lhs(F)\n", "#L = rhs(F)\n", "\n", "# Time-stepping\n", "file = File(\"./Diffusion/output.pvd\", \"compressed\")\n", "\n", "u = Function(V)\n", "t = 0\n", "for n in range(num_steps):\n", "\n", " # Update current time\n", " t += dt\n", " u_D.t = t\n", "\n", " # Compute solution\n", " solve(a == L, u, bc)\n", " \n", "\n", " # Plot solution\n", " plot(u)\n", " #plt.pause(0.5)\n", "\n", " # Compute error at vertices\n", " u_e = interpolate(u_D, V)\n", " error = np.abs(u_e.compute_vertex_values(mesh) - u.compute_vertex_values(mesh)).max()\n", " print('t = %.2f: error = %.3g' % (t, error))\n", "\n", " # Update previous solution\n", " u_n.assign(u)\n", " file << (u, t)\n", "\n", "# Hold plot\n", "#interactive()" ] } ], "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.3" } }, "nbformat": 4, "nbformat_minor": 2 }