{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Exercises\n", "\n", "\n", "\n", "\n", "\n", "## Exercise 1: Stabilizing the Crank-Nicolson method by Rannacher time stepping\n", "
\n", "\n", "mathcal{I}_t is well known that the Crank-Nicolson method may give rise to\n", "non-physical oscillations in the solution of diffusion equations\n", "if the initial data exhibit jumps (see the section [diffu:pde1:analysis:CN](#diffu:pde1:analysis:CN)).\n", "Rannacher [[Rannacher_1984]](#Rannacher_1984) suggested a stabilizing technique\n", "consisting of using the Backward Euler scheme for the first two\n", "time steps with step length $\\frac{1}{2}\\Delta t$. One can generalize\n", "this idea to taking $2m$ time steps of size $\\frac{1}{2}\\Delta t$ with\n", "the Backward Euler method and then continuing with the\n", "Crank-Nicolson method, which is of second-order in time.\n", "The idea is that the high frequencies of the initial solution are\n", "quickly damped out, and the Backward Euler scheme treats these\n", "high frequencies correctly. Thereafter, the high frequency content of\n", "the solution is gone and the Crank-Nicolson method will do well.\n", "\n", "Test this idea for $m=1,2,3$ on a diffusion problem with a\n", "discontinuous initial condition. Measure the convergence rate using\n", "the solution ([diffu:analysis:pde1:step:erf:sol](#diffu:analysis:pde1:step:erf:sol)) with the boundary\n", "conditions\n", "([diffu:analysis:pde1:p1:erf:uL](#diffu:analysis:pde1:p1:erf:uL))-([diffu:analysis:pde1:p1:erf:uR](#diffu:analysis:pde1:p1:erf:uR))\n", "for $t$ values such that the conditions are in the vicinity of $\\pm 1$.\n", "For example, $t< 5a 1.6\\cdot 10^{-2}$ makes the solution diffusion from\n", "a step to almost a straight line. The\n", "program `diffu_erf_sol.py` shows how to compute the analytical\n", "solution.\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Project 2: Energy estimates for diffusion problems\n", "
\n", "\n", "\n", "This project concerns so-called *energy estimates* for diffusion problems\n", "that can be used for qualitative analytical insight and for\n", "verification of implementations.\n", "\n", "\n", "**a)**\n", "We start with a 1D homogeneous diffusion equation with zero Dirichlet\n", "conditions:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "u_t = \\alpha u_xx, x\\in \\Omega =(0,L),\\ t\\in (0,T],\n", "\\label{diffu:exer:estimates:p1:pdf} \\tag{1} \n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "u(0,t) = u(L,t) = 0, t\\in (0,T],\n", "\\label{diffu:exer:estimates:p1:bc} \\tag{2}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "u(x,0) = I(x), x\\in [0,L]\n", "\\label{diffu:exer:estimates:p1:ic} \\tag{3}\n", "\\thinspace .\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The energy estimate for this problem reads" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "||u||_{L^2} \\leq ||I||_{L^2},\n", "\\label{diffu:exer:estimates:p1:result} \\tag{4}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where the $||\\cdot ||_{L^2}$ norm is defined by" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "||g||_{L^2} = \\sqrt{\\int_0^L g^2dx}\\thinspace .\n", "\\label{diffu:exer:estimates:L2} \\tag{5}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The quantify $||u||_{L^2}$ or $\\frac{1}{2} ||u||_{L^2}$ is known\n", "as the *energy* of the solution, although it is not the physical\n", "energy of the system. A mathematical tradition has introduced the\n", "notion *energy* in this context.\n", "\n", "The estimate ([4](#diffu:exer:estimates:p1:result)) says that the\n", "\"size of $u$\" never exceeds that of the initial condition,\n", "or more precisely, it says that the area under the $u$ curve decreases\n", "with time.\n", "\n", "To show ([4](#diffu:exer:estimates:p1:result)), multiply the PDE\n", "by $u$ and integrate from $0$ to $L$. Use that $uu_t$ can be\n", "expressed as the time derivative of $u^2$ and that $u_xxu$ can\n", "integrated by parts to form an integrand $u_x^2$. Show that\n", "the time derivative of $||u||_{L^2}^2$ must be less than or equal\n", "to zero. Integrate this expression and derive\n", "([4](#diffu:exer:estimates:p1:result)).\n", "\n", "\n", "\n", "**b)**\n", "Now we address a slightly different problem," ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "u_t = \\alpha u_xx + f(x,t), x\\in \\Omega =(0,L),\\ t\\in (0,T],\n", "\\label{diffu:exer:estimates:p2:pdf} \\tag{6} \n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "u(0,t) = u(L,t) = 0, t\\in (0,T],\n", "\\label{diffu:exer:estimates:p2:bc} \\tag{7}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "u(x,0) = 0, x\\in [0,L]\n", "\\label{diffu:exer:estimates:p2:ic} \\tag{8}\n", "\\thinspace .\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The associated energy estimate is" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "||u||_{L^2} \\leq ||f||_{L^2}\\thinspace .\n", "\\label{diffu:exer:estimates:p2:result} \\tag{9}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(This result is more difficult to derive.)\n", "\n", "Now consider the compound problem with an initial condition $I(x)$ and\n", "a right-hand side $f(x,t)$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "u_t = \\alpha u_xx + f(x,t), x\\in \\Omega =(0,L),\\ t\\in (0,T],\n", "\\label{diffu:exer:estimates:p3:pdf} \\tag{10} \n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "u(0,t) = u(L,t) = 0, t\\in (0,T],\n", "\\label{diffu:exer:estimates:p3:bc} \\tag{11}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "u(x,0) = I(x), x\\in [0,L]\n", "\\label{diffu:exer:estimates:p3:ic} \\tag{12}\n", "\\thinspace .\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Show that if $w_1$ fulfills\n", "([1](#diffu:exer:estimates:p1:pdf))-([3](#diffu:exer:estimates:p1:ic))\n", "and $w_2$ fulfills\n", "([6](#diffu:exer:estimates:p2:pdf))-([8](#diffu:exer:estimates:p2:ic)),\n", "then $u=w_1 + w_2$ is the solution of\n", "([10](#diffu:exer:estimates:p3:pdf))-([12](#diffu:exer:estimates:p3:ic)).\n", "Using the triangle inequality for norms," ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "||a + b|| \\leq ||a|| + ||b||,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "show that the energy estimate for\n", "([10](#diffu:exer:estimates:p3:pdf))-([12](#diffu:exer:estimates:p3:ic))\n", "becomes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "||u||_{L^2} \\leq ||I||_{L^2} + ||f||_{L^2}\\thinspace .\n", "\\label{diffu:exer:estimates:p3:result} \\tag{13}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**c)**\n", "One application of ([13](#diffu:exer:estimates:p3:result)) is to prove uniqueness\n", "of the solution.\n", "Suppose $u_1$ and $u_2$ both fulfill\n", "([10](#diffu:exer:estimates:p3:pdf))-([12](#diffu:exer:estimates:p3:ic)).\n", "Show that $u=u_1 - u_2$ then fulfills\n", "([10](#diffu:exer:estimates:p3:pdf))-([12](#diffu:exer:estimates:p3:ic))\n", "with $f=0$ and $I=0$. Use ([13](#diffu:exer:estimates:p3:result))\n", "to deduce that the energy must be zero for all times and therefore\n", "that $u_1=u_2$, which proves that the solution is unique.\n", "\n", "**d)**\n", "Generalize ([13](#diffu:exer:estimates:p3:result)) to a 2D/3D\n", "diffusion equation $u_t = \\nabla\\cdot (\\alpha \\nabla u)$ for $x\\in\\Omega$.\n", "\n", "\n", "\n", "**Hint.**\n", "Use integration by parts in multi dimensions:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\int_\\Omega u \\nabla\\cdot (\\alpha\\nabla u)\\dx =\n", "- \\int_\\Omega \\alpha \\nabla u\\cdot\\nabla u\\dx\n", "+ \\int_{\\partial\\Omega} u \\alpha\\frac{\\partial u}{\\partial n},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $\\frac{\\partial u}{\\partial n} = \\boldsymbol{n}\\cdot\\nabla u$,\n", "$\\boldsymbol{n}$ being the outward unit normal to the boundary $\\partial\\Omega$\n", "of the domain $\\Omega$.\n", "\n", "\n", "\n", "**e)**\n", "Now we also consider the multi-dimensional PDE $u_t =\n", "\\nabla\\cdot (\\alpha \\nabla u)$. Integrate both sides over $\\Omega$\n", "and use Gauss' divergence theorem, $\\int_\\Omega \\nabla\\cdot\\boldsymbol{q}\\dx\n", "= \\int_{\\partial\\Omega}\\boldsymbol{q}\\cdot\\boldsymbol{n}\\ds$ for a vector field\n", "$\\boldsymbol{q}$. Show that if we have homogeneous Neumann conditions\n", "on the boundary, $\\partial u/\\partial n=0$, area under the\n", "$u$ surface remains constant in time and" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\int_{\\Omega} u\\dx = \\int_{\\Omega} I\\dx\n", "\\thinspace .\n", "\\label{diffu:exer:estimates:p4:result} \\tag{14}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**f)**\n", "Establish a code in 1D, 2D, or 3D that can solve a diffusion equation with a\n", "source term $f$, initial condition $I$, and zero Dirichlet or\n", "Neumann conditions on the whole boundary.\n", "\n", "We can use ([13](#diffu:exer:estimates:p3:result))\n", "and ([14](#diffu:exer:estimates:p4:result)) as a partial verification\n", "of the code. Choose some functions $f$ and $I$ and\n", "check that ([13](#diffu:exer:estimates:p3:result)) is obeyed at any\n", "time when zero Dirichlet conditions are used.\n", "mathcal{I}_terate over the same $I$ functions and check that\n", "([14](#diffu:exer:estimates:p4:result)) is fulfilled\n", "when using zero Neumann conditions.\n", "\n", "**g)**\n", "Make a list of some possible bugs in the code, such as indexing errors\n", "in arrays, failure to set the correct boundary conditions,\n", "evaluation of a term at a wrong time level, and similar.\n", "For each of the bugs, see if the verification tests from the previous\n", "subexercise pass or fail. This investigation shows how strong\n", "the energy estimates and the estimate ([14](#diffu:exer:estimates:p4:result))\n", "are for pointing out errors in the implementation.\n", "\n", "Filename: `diffu_energy`.\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Exercise 3: Splitting methods and preconditioning\n", "
\n", "\n", "\n", "In the section [diffu:2D:direct_vs_iter](#diffu:2D:direct_vs_iter), we outlined a class of\n", "iterative methods for $Au=b$ based on splitting $A$ into $A=M-N$\n", "and introducing the iteration" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "Mu^{k} = Nu^k + b\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The very simplest splitting is $M=I$, where $I$ is the identity\n", "matrix. Show that this choice corresponds to the iteration" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "u^k = u^{k-1} + r^{k-1},\\quad r^{k-1} = b - Au^{k-1},\n", "\\label{diffu:exer:splitting_prec:simplest} \\tag{15}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $r^{k-1}$ is the residual in the linear system in iteration\n", "$k-1$. The formula ([15](#diffu:exer:splitting_prec:simplest)) is known\n", "as Richardson's iteration.\n", "Show that if we apply the simple iteration method\n", "([15](#diffu:exer:splitting_prec:simplest)) to the *preconditioned*\n", "system $M^{-1}Au=M^{-1}b$, we arrive at the Jacobi method by choosing\n", "$M=D$ (the diagonal of $A$) as preconditioner and the SOR method by\n", "choosing $M=\\omega^{-1}D + L$ ($L$ being the lower triangular part of\n", "$A$). This equivalence shows that we can apply one iteration of the\n", "Jacobi or SOR method as preconditioner.\n", "\n", "\n", "\n", "**Solution.**\n", "Inserting $M=I$ and $N=I-A$ in the iterative method leads to" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u^{k} = (I-A)u^{k-1} + b = u^{k-1} + (b - Au^{k-1}),\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "which is ([15](#diffu:exer:splitting_prec:simplest)).\n", "Replacing $A$ by $M^{-1}A$ and $b$ by $M^{-1}b$ in this equation\n", "gives" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u^k = u^{k-1} + M^{-1}r^{k-1},\\quad r^{k-1}=b-Au^{k-1},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "which we after multiplication by $M$ and reordering can write\n", "as" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "Mu^k = (M-A)u^{k-1} + b = Nu^{k-1} + b,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "which is the standard form for the Jacobi and SOR methods. Choosing $M=D$\n", "gives Jacobi and $M=\\omega^{-1}D+L$ gives SOR. We have shown that we may\n", "view $M$ as a preconditioner of a simplest possible iteration method.\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Problem 4: Oscillating surface temperature of the earth\n", "
\n", "\n", "Consider a day-and-night or seasonal variation in temperature at\n", "the surface of the earth. How deep down in the ground will the\n", "surface oscillations reach? For simplicity, we model only the\n", "vertical variation along a coordinate $x$, where $x=0$ at the\n", "surface, and $x$ increases as we go down in the ground.\n", "The temperature is governed by the heat equation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\varrho c_v\\frac{\\partial T}{\\partial t} = \\nabla\\cdot(k\\nabla T),\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "in some spatial domain $x\\in [0,L]$, where $L$ is chosen large enough such\n", "that we can assume that $T$ is approximately constant, independent of the surface\n", "oscillations, for $x>L$. The parameters $\\varrho$, $c_v$, and $k$ are the\n", "density, the specific heat capacity at constant volume, and the\n", "heat conduction coefficient, respectively.\n", "\n", "\n", "**a)**\n", "Derive the mathematical model for computing $T(x,t)$.\n", "Assume the surface oscillations to be sinusoidal around some mean\n", "temperature $T_m$. Let $T=T_m$ initially. At $x=L$, assume $T\\approx T_m$.\n", "\n", "\n", "\n", "**Solution.**\n", "The surface temperature is set as" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "T(0,t) = T_m + A\\sin(\\omega t)\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With only one \"active\" spatial coordinate we get the initial-boundary\n", "value problem" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{alignat*}{2}\n", "\\varrho c_v \\frac{\\partial T}{\\partial t} &= \\frac{\\partial}{\\partial x}\n", "\\left(k(x)\\frac{\\partial T}{\\partial x}\\right), & x\\in (0,L),\\ t\\in (0,T],\\\\\n", "T(x,0)&= T_m, & x\\in [0,L],\\\\\n", "T(0,t)&= T_m + A\\sin(\\omega t), & t\\in (0,T],\\\\\n", "T(L,t) &= T_m, & t\\in (0,T].\n", "\\end{alignat*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "**b)**\n", "Scale the model in a) assuming $k$ is constant. Use a time scale\n", "$t_c = \\omega^{-1}$ and a length scale $x_c = \\sqrt{2\\dfc/\\omega}$,\n", "where $\\dfc = k/(\\varrho c_v)$. The primary unknown can be scaled\n", "as $\\frac{T-T_m}{2A}$.\n", "\n", "Show that the scaled PDE is" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{\\partial u}{\\partial \\bar t} =\n", "\\frac{1}{2}\\frac{\\partial^2 u}{\\partial x^2},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "with initial condition $u(\\bar x,0) = 0$,\n", "left boundary condition\n", "$u(0,\\bar t) = \\sin(\\bar t)$,\n", "and right boundary condition\n", "$u(\\bar L,\\bar t) = 0$. The bar indicates a dimensionless quantity.\n", "\n", "Show that $u(\\bar x, \\bar t)=e^{-\\bar x}\\sin (\\bar x - \\bar t)$ is a\n", "solution that fulfills the PDE and the boundary condition at $\\bar x\n", "=0$ (this is the solution we will experience as $\\bar\n", "t\\rightarrow\\infty$ and $L\\rightarrow\\infty$). Conclude that an\n", "appropriate domain for $x$ is $[0,4]$ if a damping $e^{-4}\\approx\n", "0.18$ is appropriate for implementing $\\bar u\\approx\\hbox{const}$;\n", "increasing to $[0,6]$ damps $\\bar u$ to 0.0025.\n", "\n", "\n", "\n", "**Solution.**\n", "Chapter 3.2.4 in the book [[Langtangen_scaling]](#Langtangen_scaling) describes the\n", "scaling of this problem in detail.\n", "Inserting dimensionless variables $\\bar t = \\omega t$, $\\bar x =\n", "\\sqrt{\\omega/(2\\dfc)} x$, and" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u = \\frac{T-T_m}{2A},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "leads to" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{alignat*}{2}\n", "\\frac{\\partial u}{\\partial \\bar t} &=\n", "\\frac{1}{2}\\frac{\\partial^2 u}{\\partial x^2},\n", "\\quad & \\bar x\\in (0,\\bar L),\\ \\bar t\\in (0,\\bar T],\n", "\\\\\n", "u(\\bar x,0) &= 0,\n", "\\quad &\\bar x\\in [0,1],\n", "\\\\\n", "u(0,\\bar t) & = \\sin(\\bar t),\n", "\\quad &\\bar t\\in (0,\\bar T],\n", "\\\\\n", "u(\\bar L,\\bar t) & = 0,\n", "\\quad &\\bar t\\in (0,\\bar T].\n", "\\end{alignat*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The domain lengths $\\bar L$ and $\\bar T$ follows from straightforward\n", "scaling of $L$ and $T$.\n", "\n", "Inserting $u(\\bar x, \\bar t)=e^{-\\bar x}\\sin (\\bar t - \\bar x)$ in the\n", "PDE shows that this is a solution. mathcal{I}_t also obeys\n", "the boundary condition $\\bar u(0,\\bar t)=sin(\\bar t)$. As\n", "$\\bar t\\rightarrow\\infty$, the initial condition has no longer impact\n", "on the solution and is \"forgotten\" and of no interest.\n", "The boundary condition at $\\bar x=\\bar L$ is never compatible with the\n", "given solution unless $\\bar u$ is damped to zero, which happens\n", "mathematically as $\\bar L\\rightarrow\\infty$. For a numerical solution,\n", "however, we may use a small finite value such as $\\bar L=4$.\n", "\n", "\n", "\n", "**c)**\n", "Compute the scaled temperature and make animations comparing two solutions\n", "with $\\bar L=4$ and $\\bar L=8$, respectively (keep $\\Delta x$ the same).\n", "\n", "\n", "\n", "**Solution.**\n", "We can use the `viz` function in `diff1D_vc.py` to do the number\n", "crunching. Appropriate calls and visualization go here:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "\n", "import sys, os\n", "sys.path.insert(0, os.path.join(os.pardir, 'src-diffu'))\n", "from diffu1D_vc import viz\n", "\n", "sol = [] # store solutions\n", "for Nx, L in [[20, 4], [40, 8]]:\n", " dt = 0.1\n", " dx = float(L)/Nx\n", " D = dt/dx**2\n", " from math import pi, sin\n", " T = 2*pi*6\n", " from numpy import zeros\n", " a = zeros(Nx+1) + 0.5\n", " cpu, u_ = viz(\n", " I=lambda x: 0, a=a, L=L, Nx=Nx, D=D, T=T,\n", " umin=-1.1, umax=1.1, theta=0.5,\n", " u_L=lambda t: sin(t),\n", " u_R=0,\n", " animate=False, store_u=True)\n", " sol.append(u_)\n", " print('computed solution for Nx=%d in [0,%g]' % (Nx, L))\n", "\n", "print sol[0].shape\n", "print sol[1].shape\n", "import scitools.std as plt\n", "counter = 0\n", "for u0, u1 in zip(sol[0][2:], sol[1][2:]):\n", " x0 = sol[0][0]\n", " x1 = sol[1][0]\n", " plt.plot(x0, u0, 'r-', x1, u1, 'b-',\n", " legend=['short', 'long'],\n", " savefig='tmp_%04d.png' % counter,\n", " axis=[x1[0], x1[-1], -1.1, 1.1])\n", " counter += 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from IPython.display import HTML\n", "_s = \"\"\"\n", "
\n", "\n", "
\n", "

\n", "\n", "\n", "\n", "\n", "\"\"\"\n", "HTML(_s)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Problem 5: Oscillating and pulsating flow in tubes\n", "
\n", "\n", "We consider flow in a straight tube with radius $R$ and straight walls.\n", "The flow is driven by a pressure gradient $\\beta(t)$. The effect of\n", "gravity can be neglected. The mathematical problem reads" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\varrho\\frac{\\partial u}{\\partial t} =\n", "\\mu\\frac{1}{r}\\frac{\\partial}{\\partial r}\\left(\n", "r\\frac{\\partial u}{\\partial r}\\right) + \\beta(t),\\quad\n", " r\\in [0,R],\\ t\\in (0,T],\n", "\\label{_auto1} \\tag{16}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "u(r,0) = I(r),\\quad r\\in [0,R],\n", "\\label{_auto2} \\tag{17}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "u(R,t) = 0,\\quad t\\in (0,T],\n", "\\label{_auto3} \\tag{18}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "\\frac{\\partial u}{\\partial r}(0,t) = 0,\\quad t\\in (0,T].\n", "\\label{_auto4} \\tag{19}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We consider two models for $\\beta(t)$. One plain, sinusoidal oscillation:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\beta = A\\sin(\\omega t),\n", "\\label{_auto5} \\tag{20}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and one with periodic pulses," ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\beta = A\\sin^{16}(\\omega t),\n", "\\label{_auto6} \\tag{21}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that both models can be written as $\\beta = A\\sin^m(\\omega t)$, with\n", "$m=1$ and $m=16$, respectively.\n", "\n", "\n", "**a)**\n", "Scale the mathematical model, using the viscous time scale $\\varrho R^2/\\mu$.\n", "\n", "\n", "\n", "**Solution.**\n", "We can introduce" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\bar r = \\frac{r}{R}, \\quad \\bar t = \\frac{t}{\\varrho R^2/\\mu},\\quad u = \\frac{u}{u_c}\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Inserted in the PDE, we get" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{\\partial\\bar u}{\\partial\\bar t} =\n", "\\frac{1}{\\bar r}\\frac{\\partial}{\\partial\\bar r}\\left(\n", "\\bar r\\frac{\\partial\\bar u}{\\partial\\bar r}\\right) +\n", "\\frac{R^2 A}{u_c \\mu}\\sin^m (\\alpha\\bar t)\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $\\alpha$ is a dimensionless number" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\alpha = \\frac{\\omega\\varrho R^2}{\\mu} = \\frac{\\varrho R^2/\\mu}{1/\\omega},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "reflecting the ratio of the viscous diffusion time scale and the\n", "time scale of the oscillating pressure gradient.\n", "We may choose $u_c$ such that the coefficient in the pressure gradient\n", "term equals unity:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u_c = \\frac{R^2 A}{\\mu}\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The governing PDE, dropping the bars, then reads" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{\\partial u}{\\partial t} =\n", "\\frac{1}{r}\\frac{\\partial}{\\partial r}\\left(\n", "r\\frac{\\partial u}{\\partial r}\\right) +\n", "\\sin^m (\\alpha\\bar t),\\quad r\\in (0,1),\\ t\\in (0,T]\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "**b)**\n", "Implement the scaled model from a), using the unifying $\\theta$ scheme\n", "in time and centered differences in space.\n", "\n", "\n", "\n", "**Solution.**\n", "We need to take into account extensions below: a coefficient in front of\n", "the viscous term, and an extra source term.\n", "\n", "A preliminary and unfinished code:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\"\"\"\n", "Solve the diffusion equation for axi-symmetric case:\n", "\n", " u_t = 1/r * (r*a(r)*u_r)_r + f(r,t)\n", "\n", "on (0,R) with boundary conditions u(0,t)_r = 0 and u(R,t) = 0,\n", "for t in (0,T]. Initial condition: u(r,0) = I(r). \n", "Pressure gradient f.\n", "\n", "The following naming convention of variables are used.\n", "\n", "===== ==========================================================\n", "Name Description\n", "===== ==========================================================\n", "Nx The total number of mesh cells; mesh points are numbered\n", " from 0 to Nx.\n", "T The stop time for the simulation.\n", "I Initial condition (Python function of x).\n", "a Variable coefficient (constant).\n", "R Length of the domain ([0,R]).\n", "r Mesh points in space.\n", "t Mesh points in time.\n", "n Index counter in time.\n", "u Unknown at current/new time level.\n", "u_1 u at the previous time level.\n", "dr Constant mesh spacing in r.\n", "dt Constant mesh spacing in t.\n", "===== ==========================================================\n", "\n", "``user_action`` is a function of ``(u, r, t, n)``, ``u[i]`` is the\n", "solution at spatial mesh point ``r[i]`` at time ``t[n]``, where the\n", "calling code can add visualization, error computations, data analysis,\n", "store solutions, etc.\n", "\"\"\"\n", "\n", "import scipy.sparse\n", "import scipy.sparse.linalg\n", "from numpy import linspace, zeros, random, array, ones, sum, log, sqrt\n", "import time, sys\n", "import sympy as sym \n", "\n", "\n", "def solver_theta(I, a, R, Nr, D, T, theta=0.5, u_L=None, u_R=0,\n", " user_action=None, f=0):\n", " \"\"\"\n", " The array a has length Nr+1 and holds the values of\n", " a(x) at the mesh points.\n", "\n", " Method: (implicit) theta-rule in time.\n", "\n", " Nr is the total number of mesh cells; mesh points are numbered\n", " from 0 to Nr.\n", " D = dt/dr**2 and implicitly specifies the time step.\n", " T is the stop time for the simulation.\n", " I is a function of r.\n", " u_L = None implies du/dr = 0, i.e. a symmetry condition \n", " f(r,t) is pressure gradient with radius.\n", "\n", " user_action is a function of (u, x, t, n) where the calling code\n", " can add visualization, error computations, data analysis,\n", " store solutions, etc.\n", " \n", " r*alpha is needed midway between spatial mesh points, - use\n", " arithmetic mean of successive mesh values (i.e. of r_i*alpha_i)\n", " \"\"\"\n", " import time\n", " t0 = time.perf_counter()\n", "\n", " r = linspace(0, R, Nr+1) # mesh points in space\n", " dr = r[1] - r[0]\n", " dt = D*dr**2 \n", " Nt = int(round(T/float(dt)))\n", " t = linspace(0, T, Nt+1) # mesh points in time\n", "\n", " if isinstance(u_L, (float,int)):\n", " u_L_ = float(u_L) # must take copy of u_L number\n", " u_L = lambda t: u_L_\n", " if isinstance(u_R, (float,int)):\n", " u_R_ = float(u_R) # must take copy of u_R number\n", " u_R = lambda t: u_R_\n", " if isinstance(f, (float,int)):\n", " f_ = float(f) # must take copy of f number\n", " f = lambda r, t: f_\n", "\n", " ra = r*a # help array in scheme\n", "\n", " inv_r = zeros(len(r)-2) # needed for inner mesh points\n", " inv_r = 1.0/r[1:-1]\n", "\n", " u = zeros(Nr+1) # solution array at t[n+1]\n", " u_1 = zeros(Nr+1) # solution at t[n]\n", "\n", " Dl = 0.5*D*theta\n", " Dr = 0.5*D*(1-theta)\n", "\n", " # Representation of sparse matrix and right-hand side\n", " diagonal = zeros(Nr+1)\n", " lower = zeros(Nr)\n", " upper = zeros(Nr)\n", " b = zeros(Nr+1)\n", "\n", " # Precompute sparse matrix (scipy format)\n", " diagonal[1:-1] = 1 + Dl*(ra[2:] + 2*ra[1:-1] + ra[:-2])*inv_r\n", " lower[:-1] = -Dl*(ra[1:-1] + ra[:-2])*inv_r\n", " upper[1:] = -Dl*(ra[2:] + ra[1:-1])*inv_r\n", " # Insert boundary conditions\n", " if u_L == None: # symmetry axis, du/dr = 0\n", " diagonal[0] = 1 + 8*a[0]*Dl\n", " upper[0] = -8*a[0]*Dl\n", " else:\n", " diagonal[0] = 1\n", " upper[0] = 0\n", " diagonal[Nr] = 1\n", " lower[-1] = 0\n", "\n", " A = scipy.sparse.diags(\n", " diagonals=[diagonal, lower, upper],\n", " offsets=[0, -1, 1],\n", " shape=(Nr+1, Nr+1),\n", " format='csr')\n", " #print A.todense()\n", "\n", " # Set initial condition\n", " for i in range(0,Nr+1):\n", " u_1[i] = I(r[i])\n", "\n", " if user_action is not None:\n", " user_action(u_1, r, t, 0)\n", "\n", " # Time loop\n", " for n in range(0, Nt):\n", " b[1:-1] = u_1[1:-1] + Dr*(\n", " (ra[2:] + ra[1:-1])*(u_1[2:] - u_1[1:-1]) -\n", " (ra[1:-1] + ra[0:-2])*(u_1[1:-1] - u_1[:-2]))*inv_r + \\\n", " dt*theta*f(r[1:-1], t[n+1]) + \\\n", " dt*(1-theta)*f(r[1:-1], t[n])\n", " \n", " # Boundary conditions\n", " if u_L == None: # symmetry axis, du/dr = 0\n", " b[0] = u_1[0] + 8*a[0]*Dr*(u_1[1] - u_1[0]) + \\\n", " dt*theta*f(0, (n+1)*dt) + \\\n", " dt*(1 - theta)*f(0, n*dt)\n", " else: \n", " b[0] = u_L(t[n+1]) \n", " b[-1] = u_R(t[n+1])\n", " #print b \n", " \n", " # Solve\n", " u[:] = scipy.sparse.linalg.spsolve(A, b)\n", " \n", " if user_action is not None:\n", " user_action(u, r, t, n+1)\n", "\n", " # Switch variables before next step\n", " u_1, u = u, u_1\n", "\n", " t1 = time.perf_counter()\n", " # return u_1, since u and u_1 are switched\n", " return u_1, t, t1-t0\n", "\n", "def compute_rates(h_values, E_values):\n", " m = len(h_values)\n", " q = [log(E_values[i+1]/E_values[i])/\n", " log(h_values[i+1]/h_values[i])\n", " for i in range(0, m-1, 1)]\n", " q = [round(q_, 2) for q_ in q]\n", " return q\n", "\n", "def make_a(alpha, r):\n", " \"\"\"\n", " alpha is a func, generally of r, - but may be constant.\n", " Note: when solution is to be axi-symmetric, alpha\n", " must be so too.\n", " \"\"\"\n", " a = alpha(r)*ones(len(r))\n", " return a\n", "\n", "def tests_with_alpha_and_u_exact():\n", " '''\n", " Test solver performance when alpha is either const or \n", " a fu of r, combined with a manufactured sol u_exact \n", " that is either a fu of r only, or a fu of both r and t.\n", " Note: alpha and u_e are defined as symb expr here, since \n", " test_solver_symmetric needs to automatically generate \n", " the source term f. After that, test_solver_symmetric\n", " redefines alpha, u_e and f as num functions.\n", " '''\n", " R, r, t = sym.symbols('R r t')\n", "\n", " # alpha const ...\n", " \n", " # ue = const\n", " print('Testing with alpha = 1.5 and u_e = R**2 - r**2...')\n", " test_solver_symmetric(alpha=1.5, u_exact=R**2 - r**2)\n", " \n", " # ue = ue(t)\n", " print('Testing with alpha = 1.5 and u_e = 5*t*(R**2 - r**2)...')\n", " test_solver_symmetric(alpha=1.5, u_exact=5*t*(R**2 - r**2))\n", " \n", " # alpha function of r ...\n", " \n", " # ue = const \n", " print('Testing with alpha = 1 + r**2 and u_e = R**2 - r**2...')\n", " test_solver_symmetric(alpha=1+r**2, u_exact=R**2 - r**2)\n", " \n", " # ue = ue(t)\n", " print('Testing with alpha = 1+r**2 and u_e = 5*t*(R**2 - r**2)...')\n", " test_solver_symmetric(alpha=1+r**2, u_exact=5*t*(R**2 - r**2))\n", "\n", "\n", "\n", "def test_solver_symmetric(alpha, u_exact):\n", " '''\n", " Test solver performance for manufactured solution\n", " given in the function u_exact. Parameter alpha is \n", " either a const or a function of r. In the latter \n", " case, an \"exact\" sol can not be achieved, so then\n", " testing switches to conv. rates.\n", " R is tube radius and T is duration of simulation.\n", " alpha constant:\n", " Compares the manufactured solution with the \n", " solution from the solver at each time step. \n", " alpha function of r:\n", " convergence rates are tested (using the sol\n", " at the final point in time only).\n", " ''' \n", " \n", " def compare(u, r, t, n): # user_action function\n", " \"\"\"Compare exact and computed solution.\"\"\"\n", " u_e = u_exact(r, t[n])\n", " diff = abs(u_e - u).max()\n", " #print diff\n", " tol = 1E-12\n", " assert diff < tol, 'max diff: %g' % diff\n", "\n", " def pde_source_term(a, u):\n", " '''Return the terms in the PDE that the source term\n", " must balance, here du/dt - (1/r) * d/dr(r*a*du/dr).\n", " a, i.e. alpha, is either const or a fu of r.\n", " u is a symbolic Python function of r and t.'''\n", " \n", " return sym.diff(u, t) - \\\n", " (1.0/r)*sym.diff(r*a*sym.diff(u, r), r)\n", " \n", " R, r, t = sym.symbols('R r t')\n", "\n", " # fit source term\n", " f = sym.simplify(pde_source_term(alpha, u_exact)) \n", "\n", " R = 1.0 # radius of tube\n", " T = 2.0 # duration of simulation \n", " \n", " if sym.diff(alpha, r) == 0: \n", " alpha_is_const = True\n", " else:\n", " alpha_is_const = False \n", "\n", " # make alpha, f and u_exact numerical functions\n", " alpha = sym.lambdify([r], alpha, modules='numpy') \n", " f = sym.lambdify([r, t], f.subs('R', R), modules='numpy') \n", " u_exact = sym.lambdify(\n", " [r, t], u_exact.subs('R', R), modules='numpy') \n", "\n", " I = lambda r: u_exact(r, 0)\n", "\n", " # some help variables\n", " FE = 0 # Forward Euler method\n", " BE = 1 # Backward Euler method\n", " CN = 0.5 # Crank-Nicolson method\n", "\n", " # test all three schemes \n", " for theta in (FE, BE, CN):\n", " print('theta: ', theta)\n", " E_values = []\n", " dt_values = []\n", " for Nr in (2, 4, 8, 16, 32, 64):\n", " print('Nr:', Nr)\n", " r = linspace(0, R, Nr+1) # mesh points in space\n", " dr = r[1] - r[0]\n", " a_values = make_a(alpha, r) \n", " if theta == CN:\n", " dt = dr\n", " else: # either FE or BE\n", " # use most conservative dt as decided by FE\n", " K = 1.0/(4*a_values.max()) \n", " dt = K*dr**2 \n", " D = dt/dr**2\n", "\n", " if alpha_is_const: \n", " u, t, cpu = solver_theta(\n", " I, a_values, R, Nr, D, T, \n", " theta, u_L=None, u_R=0,\n", " user_action=compare, f=f) \n", " else: # alpha depends on r\n", " u, t, cpu = solver_theta(\n", " I, a_values, R, Nr, D, T, \n", " theta, u_L=None, u_R=0,\n", " user_action=None, f=f) \n", " \n", " # compute L2 error at t = T\n", " u_e = u_exact(r, t[-1])\n", " e = u_e - u\n", " E = sqrt(dr*sum(e**2))\n", " E_values.append(E)\n", " dt_values.append(dt)\n", " \n", " if alpha_is_const is False: \n", " q = compute_rates(dt_values, E_values) \n", " print('theta=%g, q: %s' % (theta, q))\n", " expected_rate = 2 if theta == CN else 1\n", " tol = 0.1\n", " diff = abs(expected_rate - q[-1])\n", " print('diff:', diff)\n", " assert diff < tol\n", " \n", " \n", "if __name__ == '__main__':\n", " tests_with_alpha_and_u_exact() \n", " print('This is just a start. More remaining for this Exerc.')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "**c)**\n", "Verify the implementation in b) using a manufactured solution that is\n", "quadratic in $r$ and linear in $t$. Make a corresponding test function.\n", "\n", "\n", "\n", "**Hint.**\n", "You need to include an extra source term\n", "in the equation to allow for such tests. Let the spatial variation be\n", "$1-r^2$ such that the boundary condition is fulfilled.\n", "\n", "\n", "\n", "**d)**\n", "Make animations for $m=1,16$ and $\\alpha=1,0.1$. Choose $T$ such that\n", "the motion has reached a steady state (non-visible changes from period to\n", "period in $u$).\n", "\n", "**e)**\n", "For $\\alpha\\gg 1$, the scaling in a) is not good, because the\n", "characteristic time for changes (due to the pressure) is much smaller\n", "than the viscous diffusion time scale ($\\alpha$ becomes large).\n", "We should in this case base\n", "the short time scale on $1/\\omega$. Scale the model again, and\n", "make an animation for $m=1,16$ and $\\alpha = 10$.\n", "\n", "\n", "\n", "**Solution.**\n", "Now the governing PDE becomes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{\\partial u}{\\partial t} =\n", "\\alpha^{-1}\\frac{1}{r}\\frac{\\partial}{\\partial r}\\left(\n", "r\\frac{\\partial u}{\\partial r}\\right) +\n", "\\sin^m t,\\quad r\\in (0,1),\\ t\\in (0,T]\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this case," ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u_c = \\frac{A}{\\varrho\\omega}\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that for $\\alpha\\gg 1$, we can neglect the viscous term, and we\n", "basically have a balance between the acceleration and the driving pressure\n", "gradient:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{\\partial u}{\\partial t} = \\sin^m t\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[hpl 1: This may be a great challenge numerically, since we have a plug\n", "independent of r that oscillates back and forth. CN is probably very\n", "unstable. Can make a point out of this. Try $\\alpha=1$ and increase\n", "gently.]\n", "\n", "\n", "\n", "Filename: `axisymm_flow`.\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Problem 6: Scaling a welding problem\n", "
\n", "\n", "Welding equipment makes a very localized heat source that moves in\n", "time. We shall investigate the heating due to welding and choose, for\n", "maximum simplicity, a one-dimensional heat equation with a fixed\n", "temperature at the ends, and we neglect melting. We shall scale the\n", "problem, and besides solving such a problem numerically, the aim is to\n", "investigate the appropriateness of alternative scalings.\n", "\n", "The governing PDE problem reads" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{alignat*}{2}\n", "\\varrho c\\frac{\\partial u}{\\partial t} &= k\\frac{\\partial^2 u}{\\partial x^2}\n", "+ f, & x\\in (0,L),\\ t\\in (0,T),\\\\\n", "u(x,0) &= U_s, & x\\in [0,L],\\\\\n", "u(0,t) = u(L,t) &= 0, & t\\in (0,T].\n", "\\end{alignat*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here, $u$ is the temperature, $\\varrho$ the density of the material,\n", "$c$ a heat capacity, $k$ the heat conduction coefficient, $f$ is\n", "the heat source from the welding equipment, and $U_s$ is the\n", "initial constant (room) temperature in the material.\n", "\n", "A possible model for the heat source is a moving Gaussian function:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "f = A\\exp{\\left(-\\frac{1}{2}\\left(\\frac{x-vt}{\\sigma}\\right)^2\\right)},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $A$ is the strength, $\\sigma$ is a parameter governing how\n", "peak-shaped (or localized in space) the heat source is, and\n", "$v$ is the velocity (in positive $x$ direction) of the source.\n", "\n", "\n", "**a)**\n", "Let $x_c$, $t_c$, $u_c$, and $f_c$ be scales, i.e., characteristic\n", "sizes, of $x$, $t$, $u$, and $f$, respectively. The natural choice of\n", "$x_c$ and $f_c$ is $L$ and $A$, since these make the scaled $x$ and\n", "$f$ in the interval $[0,1]$. If each of the three terms in the PDE\n", "are equally important, we can find $t_c$ and $u_c$ by demanding that\n", "the coefficients in the scaled PDE are all equal to unity. Perform\n", "this scaling. Use scaled quantities in the arguments for the\n", "exponential function in $f$ too and show that" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\bar f= e^{-\\frac{1}{2}\\beta^2(\\bar x -\\gamma \\bar t)^2},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $\\beta$ and $\\gamma$ are dimensionless numbers. Give an\n", "interpretation of $\\beta$ and $\\gamma$.\n", "\n", "\n", "\n", "**Solution.**\n", "We introduce" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\bar x=\\frac{x}{L},\\quad \\bar t = \\frac{t}{t_c},\\quad \\bar u = \\frac{u-U_s}{u_c},\n", "\\quad \\bar f=\\frac{f}{A}\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Inserted in the PDE and dividing by $\\varrho c u_c/t_c$ such that the\n", "coefficient in front of $\\partial\\bar u/\\partial\\bar t$ becomes unity,\n", "and thereby all terms become dimensionless, we get" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{\\partial\\bar u}{\\partial\\bar t} =\n", "\\frac{k t_c}{\\varrho c L^2}\\frac{\\partial^2\\bar u}{\\partial\\bar x^2}\n", "+ \\frac{A t_c}{\\varrho c u_c}\\bar f\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Demanding that all three terms are equally important, it follows that" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{k t_c}{\\varrho c L^2} = 1,\\quad \\frac{A t_c}{\\varrho c u_c}=1\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These constraints imply the *diffusion time scale*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "t_c = \\frac{\\varrho cL^2}{k},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and a scale for $u_c$," ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u_c = \\frac{AL^2}{k}\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The scaled PDE reads" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{\\partial\\bar u}{\\partial\\bar t} =\n", "\\frac{\\partial^2\\bar u}{\\partial\\bar x^2}\n", "+ \\bar f\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Scaling $f$ results in" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "\\bar f &= \\exp{\\left(-\\frac{1}{2}\\left(\\frac{x-vt}{\\sigma}\\right)^2\\right)}\\\\\n", "&= \\exp{\\left(-\\frac{1}{2}\\frac{L^2}{\\sigma^2}\n", "\\left(\\bar x- \\frac{vt_c}{L}t\\right)^2\\right)}\\\\\n", "&= \\exp{\\left(-\\frac{1}{2}\\beta^2\\left(\\bar x-\\gamma \\bar t\\right)^2\\right)},\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $\\beta$ and $\\gamma$ are dimensionless numbers:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\beta = \\frac{L}{\\sigma},\\quad\n", "\\gamma = \\frac{vt_c}{L} = \\frac{v\\varrho cL}{k}\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The $\\sigma$ parameter measures the width of the Gaussian peak, so\n", "$\\beta$ is the ratio of the domain and the width of the heat source (large\n", "$\\beta$ implies a very peak-formed heat source). The $\\gamma$\n", "parameter arises from $t_c/(L/v)$, which is the ratio of the diffusion\n", "time scale and the time it takes for the heat source to travel through\n", "the domain. Equivalently, we can multiply by $t_c/t_c$ to get $\\gamma\n", "= v/(t_cL)$ as the ratio between the velocity of the heat source and\n", "the diffusion velocity.\n", "\n", "\n", "\n", "**b)**\n", "Argue that for large $\\gamma$ we should base the time scale on the\n", "movement of the heat source. Show that this gives rise to the scaled\n", "PDE" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{\\partial\\bar u}{\\partial\\bar t} =\n", "\\gamma^{-1}\\frac{\\partial^2\\bar u}{\\partial\\bar x^2}\n", "+ \\bar f,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\bar f = \\exp{(-\\frac{1}{2}\\beta^2(\\bar x - \\bar t)^2)}\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Discuss when the scalings in a) and b) are appropriate.\n", "\n", "\n", "\n", "**Solution.**\n", "We perform the scaling as in a), but this time we determine $t_c$ such\n", "that the heat source moves with unit velocity. This means that" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{vt_c}{L} = 1\\quad\\Rightarrow\\quad t_c = \\frac{L}{v}\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Scaling of the PDE gives, as before," ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{\\partial\\bar u}{\\partial\\bar t} =\n", "\\frac{k t_c}{\\varrho c L^2}\\frac{\\partial^2\\bar u}{\\partial\\bar x^2}\n", "+ \\frac{A t_c}{\\varrho c u_c}\\bar f\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Inserting the expression for $t_c$, we have" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{\\partial\\bar u}{\\partial\\bar t} =\n", "\\frac{k L}{\\varrho c L^2v}\\frac{\\partial^2\\bar u}{\\partial\\bar x^2}\n", "+ \\frac{A L}{v\\varrho c u_c}\\bar f\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We recognize the first coefficient as $\\gamma^{-1}$, while $u_c$ can\n", "be determined from demanding the second coefficient to be unity:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u_c = \\frac{AL}{v\\varrho c}\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The scaled PDE is therefore" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{\\partial\\bar u}{\\partial\\bar t} =\n", "\\gamma^{-1}\\frac{\\partial^2\\bar u}{\\partial\\bar x^2}\n", "+ \\bar f\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If the heat source moves very fast, there is little time for the\n", "diffusion to transport the heat away from the source, and the heat\n", "conduction term becomes insignificant. This is reflected in the\n", "coefficient $\\gamma^{-1}$, which is small when $\\gamma$, the ratio of\n", "the heat source velocity and the diffusion velocity, is large.\n", "\n", "The scaling in a) is therefore appropriate if diffusion is a\n", "significant process, i.e., the welding equipment moves at a slow speed\n", "so heat can efficiently spread out by diffusion. For large $\\gamma$,\n", "the scaling in b) is appropriate, and $t=1$ corresponds to having the\n", "heat source traveled through the domain (with the scaling in a), the\n", "heat source will leave the domain in short time).\n", "\n", "\n", "\n", "**c)**\n", "One aim with scaling is to get a solution that lies in the interval\n", "$[-1,1]$. This is not always the case when $u_c$ is based on a scale\n", "involving a source term, as we do in a) and b). However, from the\n", "scaled PDE we realize that if we replace $\\bar f$ with $\\delta\\bar f$,\n", "where $\\delta$ is a dimensionless factor, this corresponds to\n", "replacing $u_c$ by $u_c/\\delta$. So, if we observe that $\\bar\n", "u\\sim1/\\delta$ in simulations, we can just replace $\\bar f$ by $\\delta\n", "\\bar f$ in the scaled PDE.\n", "\n", "Use this trick and implement the two scaled models. Reuse software for\n", "the diffusion equation (e.g., the `solver` function in\n", "`diffu1D_vc.py`). Make a function `run(gamma, beta=10, delta=40,\n", "scaling=1, animate=False)` that runs the model with the given\n", "$\\gamma$, $\\beta$, and $\\delta$ parameters as well as an indicator\n", "`scaling` that is 1 for the scaling in a) and 2 for the scaling in\n", "b). The last argument can be used to turn screen animations on or off.\n", "\n", "Experiments show that with $\\gamma=1$ and $\\beta=10$, $\\delta =20$\n", "is appropriate. Then $\\max |\\bar u|$ will be larger than 4 for $\\gamma\n", "=40$, but that is acceptable.\n", "\n", "Equip the `run` function with visualization, both animation of $\\bar u$\n", "and $\\bar f$, and plots with $\\bar u$ and $\\bar f$ for $t=0.2$ and $t=0.5$.\n", "\n", "\n", "\n", "**Hint.**\n", "Since the amplitudes of $\\bar u$ and $\\bar f$ differs by a factor $\\delta$,\n", "it is attractive to plot $\\bar f/\\delta$ together with $\\bar u$.\n", "\n", "\n", "\n", "\n", "\n", "**Solution.**\n", "Here is a possible `run` function:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# from .diffu1D_vc import solver\n", "import numpy as np\n", "\n", "def run(gamma, beta=10, delta=40, scaling=1, animate=False):\n", " \"\"\"Run the scaled model for welding.\"\"\"\n", " if scaling == 1:\n", " v = gamma\n", " a = 1\n", " elif scaling == 2:\n", " v = 1\n", " a = 1.0/gamma\n", "\n", " b = 0.5*beta**2\n", " L = 1.0\n", " ymin = 0\n", " # Need gloal to be able change ymax in closure process_u\n", " global ymax\n", " ymax = 1.2\n", "\n", " I = lambda x: 0\n", " f = lambda x, t: delta*np.exp(-b*(x - v*t)**2)\n", "\n", " import time\n", " import scitools.std as plt\n", " plot_arrays = []\n", "\n", " def process_u(u, x, t, n):\n", " global ymax\n", " if animate:\n", " plt.plot(x, u, 'r-',\n", " x, f(x, t[n])/delta, 'b-',\n", " axis=[0, L, ymin, ymax], title='t=%f' % t[n],\n", " xlabel='x', ylabel='u and f/%g' % delta)\n", " if t[n] == 0:\n", " time.sleep(1)\n", " plot_arrays.append(x)\n", " dt = t[1] - t[0]\n", " tol = dt/10.0\n", " if abs(t[n] - 0.2) < tol or abs(t[n] - 0.5) < tol:\n", " plot_arrays.append((u.copy(), f(x, t[n])/delta))\n", " if u.max() > ymax:\n", " ymax = u.max()\n", "\n", " Nx = 100\n", " D = 10\n", " T = 0.5\n", " u_L = u_R = 0\n", " theta = 1.0\n", " cpu = solver(\n", " I, a, f, L, Nx, D, T, theta, u_L, u_R, user_action=process_u)\n", " x = plot_arrays[0]\n", " plt.figure()\n", " for u, f in plot_arrays[1:]:\n", " plt.plot(x, u, 'r-', x, f, 'b--', axis=[x[0], x[-1], 0, ymax],\n", " xlabel='$x$', ylabel=r'$u, \\ f/%g$' % delta)\n", " plt.hold('on')\n", " plt.legend(['$u,\\\\ t=0.2$', '$f/%g,\\\\ t=0.2$' % delta,\n", " '$u,\\\\ t=0.5$', '$f/%g,\\\\ t=0.5$' % delta])\n", " filename = 'tmp1_gamma%g_s%d' % (gamma, scaling)\n", " s = 'diffusion' if scaling == 1 else 'source'\n", " plt.title(r'$\\beta = %g,\\ \\gamma = %g,\\ $' % (beta, gamma)\n", " + 'scaling=%s' % s)\n", " plt.savefig(filename + '.pdf'); plt.savefig(filename + '.png')\n", " return cpu" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that we have dropped the bar notation in the plots. mathcal{I}_t is common\n", "to drop the bars as soon as the scaled problem is established.\n", "\n", "\n", "\n", "**d)**\n", "Use the software in c) to investigate $\\gamma=0.2,1,5,40$ for the\n", "two scalings. Discuss the results.\n", "\n", "\n", "\n", "**Solution.**\n", "For these investigations, we compare the two scalings for each of\n", "the different $\\gamma$ values. An appropriate function for automating\n", "the tasks is" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def investigate():\n", " \"\"\"Do scienfic experiments with the run function above.\"\"\"\n", " # Clean up old files\n", " import glob\n", " for filename in glob.glob('tmp1_gamma*') + \\\n", " glob.glob('welding_gamma*'):\n", " os.remove(filename)\n", "\n", " gamma_values = 1, 40, 5, 0.2, 0.025\n", " for gamma in gamma_values:\n", " for scaling in 1, 2:\n", " run(gamma=gamma, beta=10, delta=20, scaling=scaling)\n", "\n", " # Combine images\n", " for gamma in gamma_values:\n", " for ext in 'pdf', 'png':\n", " cmd = 'doconce combine_images -2 '\\\n", " 'tmp1_gamma%(gamma)g_s1.%(ext)s '\\\n", " 'tmp1_gamma%(gamma)g_s2.%(ext)s '\\\n", " 'welding_gamma%(gamma)g.%(ext)s' % vars()\n", " os.system(cmd)\n", " # pdflatex doesn't like 0.2 in filenames...\n", " if '.' in str(gamma):\n", " os.rename(\n", " 'welding_gamma%(gamma)g.%(ext)s' % vars(),\n", " ('welding_gamma%(gamma)g' % vars()).replace('.', '_')\n", " + '.' + ext)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We run here a Backward Euler scheme with $N_x=100$ and quite long\n", "time steps.\n", "\n", "Running the `investigate` function, we get the following plots:\n", "\n", "\n", "\n", "\n", "

\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "

\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "

\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "

\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "

\n", "\n", "\n", "\n", "\n", "\n", "For $\\gamma\\ll 1$ as in $\\gamma = 0.025$, the heat source moves very\n", "slowly on the diffusion time scale and has hardly entered the medium,\n", "while the scaling in b) is not inappropriate, but a larger $\\delta$ is\n", "needed to bring $\\bar u$ around unity. We see that for $\\gamma=0.2$,\n", "each of the scalings work, but with the diffusion time scale, the heat\n", "source has not moved much into the domain. For $\\gamma=1$, the\n", "mathematical problems are identical and hence the plots too. For\n", "$\\gamma=5$, the time scale based on the source is clearly the best\n", "choice, and for $\\gamma=40$, only this scale is appropriate.\n", "\n", "A conclusion is that the scaling in b) works well for a range of $\\gamma$\n", "values, also in the case $\\gamma\\ll 1$.\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Filename: `welding`.\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Exercise 7: Implement a Forward Euler scheme for axi-symmetric diffusion\n", "
\n", "\n", "Based on the discussion in the section [diffu:fd2:radial](#diffu:fd2:radial), derive in detail\n", "the discrete equations for a Forward Euler in time, centered in space,\n", "finite difference method for axi-symmetric diffusion. The\n", "diffusion coefficient may be a function of the radial coordinate.\n", "At the outer boundary $r=R$, we may have either a Dirichlet or Robin\n", "condition.\n", "Implement this scheme. Construct appropriate test problems.\n", "\n", "\n", "\n", "**Solution.**\n", "We start with the equation at $r=0$. According to the section [diffu:fd2:radial](#diffu:fd2:radial),\n", "we get" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{u^{n+1}_0-u^n_0}{\\Delta t} = 4\\dfc(0)\\frac{u_1^n - u^n_0}{\\Delta r^2}\n", "+ f_0^n\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For $i>0$, we have" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "\\frac{u^{n+1}_i-u^n_i}{\\Delta t} &= \\frac{1}{r_i\\Delta r^2}(\n", "\\frac{1}{2}(r_i + r_{i+1})\\frac{1}{2}(\\dfc_i + \\dfc_{i+1})(u^n_{i+1} - u^n_i) -\\\\\n", "&\\qquad\\frac{1}{2}(r_{i-1} + r_{i})\\frac{1}{2}(\\dfc_{i-1} + \\dfc_{i})(u^n_{i} - u^n_{i-1}))\n", "+ f_i^n\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solving with respect to $u^{n+1}_i$ and introducing $D=\\Delta t/\\Delta r^2$\n", "results in" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "u^{n+1}_0 &= u^n_0 + 4D\\dfc(0)(u_1^n - u^n_0)\n", "+ f_0^n,\\\\\n", "u^{n+1}_i &= u^n_i + D\\frac{1}{r_i}(\n", "\\frac{1}{2}(r_i + r_{i+1})\\frac{1}{2}(\\dfc_i + \\dfc_{i+1})(u^n_{i+1} - u^n_i) -\\\\\n", "&\\qquad\\frac{1}{2}(r_{i-1} + r_{i})\\frac{1}{2}(\\dfc_{i-1} + \\dfc_{i})(u^n_{i} - u^n_{i-1}))\n", "+ \\Delta t f_i^n,\\\\\n", "&\\qquad i = 1,\\ldots,N_r-1,\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and $u^{n+1}_i$ at the end point $i=N_r$ is assumed known in case of\n", "a Dirichlet condition. A Robin condition" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "-\\dfc\\frac{\\partial u}{\\partial n} = h_T(u-U_s),\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "can be discretized at $i=N_r$ by" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "-\\alpha_i\\frac{u_{i+1}^n-u_{i-1}^n}{2\\Delta r} = h_T(u_i^n - U_s)\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solving with respect to the value at the fictitious point $i+1$ gives" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u_{i+1}^n = u_{i-1}^n - 2\\Delta r \\frac{h_T}{\\alpha_i}(u_i^n - U_s)\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This value is then inserted for $u_{i+1}^n$ in the discrete PDE at $i=N_r$.\n", "\n", "\n", "Filename: `FE_axisym`.\n", "\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.8.2" } }, "nbformat": 4, "nbformat_minor": 4 }