{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# 3.1 Parabolic model problem\n",
    "In this unit we want to turn to instationary problems. And we will start with a very basic setup: implicit Euler time stepping for the convection diffusion equation. After the main part of this tutorial we have supplementary material to extend the case in several regards.\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "\n",
    "We are solving the unsteady heat equation \n",
    "\n",
    "$$\\text{find } u:[0,T] \\to H_{0,D}^1 \\quad \\int_{\\Omega} \\partial_t u v + \\int_{\\Omega} \\nabla u \\nabla v + b \\cdot \\nabla u v = \\int f v  \\quad \\forall v \\in H_{0,D}^1, \\quad u(t=0) = u_0$$\n",
    "\n",
    "with a suitable advective field $b$ (the wind)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "#imports\n",
    "from ngsolve import *\n",
    "from netgen.geom2d import SplineGeometry\n",
    "from ngsolve.webgui import Draw"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "* Geometry: $(-1,1)^2$\n",
    "* Dirichlet boundaries everywhere\n",
    "* Mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "from netgen.occ import *\n",
    "from netgen.webgui import Draw as DrawGeo\n",
    "shape = Rectangle(2,2).Face().Move((-1,-1,0))\n",
    "shape.edges.Min(X).name=\"left\"\n",
    "shape.edges.Max(X).name=\"right\"\n",
    "shape.edges.Min(Y).name=\"bottom\"\n",
    "shape.edges.Max(Y).name=\"top\"\n",
    "mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.25))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = H1(mesh, order=3, dirichlet=\"bottom|right|left|top\")\n",
    "u,v = fes.TnT()\n",
    "time = 0.0\n",
    "dt = 0.001"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We define the field $b$ (the wind) as \n",
    "\n",
    "$$b(x,y) := (2y(1-x^2),-2x(1-y^2)).$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "b = CoefficientFunction((2*y*(1-x*x),-2*x*(1-y*y)))\n",
    "Draw(b,mesh,\"wind\", vectors={\"grid_size\": 32}, order=3)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* bilinear forms for \n",
    " * convection-diffusion stiffness and \n",
    " * mass matrix separately.\n",
    "* non-symmetric memory layout for the mass matrix so that a and m have the same sparsity pattern."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "a = BilinearForm(fes, symmetric=False)\n",
    "a += 0.01*grad(u)*grad(v)*dx + b*grad(u)*v*dx\n",
    "a.Assemble()\n",
    "\n",
    "m = BilinearForm(fes, symmetric=False)\n",
    "m += u*v*dx\n",
    "m.Assemble()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Implicit Euler method\n",
    "$$\n",
    "  \\underbrace{(M + \\Delta t A)}_{M^\\ast} u^{n+1} = M u^n + \\Delta tf^{n+1}\n",
    "$$\n",
    "\n",
    "or in an incremental form:\n",
    "\n",
    "$$\n",
    "  M^\\ast (u^{n+1} - u^n) = - \\Delta t A u^n + \\Delta tf^{n+1}.\n",
    "$$\n",
    "\n",
    "* Incremental form: $u^{n+1} - u^n$ has homogeneous boundary conditions (unless boundary conditions are time-dependent).\n",
    "* For the time stepping method: set up linear combinations of matrices.\n",
    "* (Only) if the sparsity pattern of the matrices agree we can copy the pattern and sum up the entries."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "First, we create a matrix of same format as m.mat and compare number of non-zero entries and sparsity pattern:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "mstar = m.mat.CreateMatrix()\n",
    "print(f\"m.mat.nze = {m.mat.nze}, a.mat.nze={a.mat.nze}, mstar.nze={mstar.nze}\")\n",
    "from helper import ShowPattern\n",
    "print(\"sparsity pattern a.mat:\")\n",
    "ShowPattern(a.mat,binarize=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": [
     "outputPrepend"
    ]
   },
   "outputs": [],
   "source": [
    "print(\"sparsity pattern mstar:\")\n",
    "ShowPattern(mstar,binarize=True)\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "To access the entries we use the vector of nonzero-entries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "print(f\"mstar.nze={mstar.nze}, len(mstar.AsVector())={len(mstar.AsVector())}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(mstar.AsVector())"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Using the vector we can build the linear combination of the a and the m matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "mstar.AsVector().data = m.mat.AsVector() + dt * a.mat.AsVector()\n",
    "# corresponds to M* = M + dt * A\n",
    "invmstar = mstar.Inverse(freedofs=fes.FreeDofs())"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We set the r.h.s. $f = exp(-6 ((x+\\frac12)^2+y^2)) - exp(-6 ((x-\\frac12)^2+y^2))$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "f = LinearForm(fes)\n",
    "gaussp = exp(-6*((x+0.5)*(x+0.5)+y*y))-exp(-6*((x-0.5)*(x-0.5)+y*y))\n",
    "Draw(gaussp,mesh,\"f\", deformation=True)\n",
    "f += gaussp*v*dx\n",
    "f.Assemble()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "and the initial data: $u_0 = (1-y^2)x$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "gfu = GridFunction(fes)\n",
    "gfu.Set((1-y*y)*x) # note that boundary conditions remain\n",
    "scene = Draw(gfu,mesh,\"u\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "we define a simple time loop including some visualization sampling:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "def TimeStepping(invmstar, initial_cond = None, t0 = 0, tend = 2, \n",
    "                 nsamples = 10):\n",
    "    if initial_cond:\n",
    "        gfu.Set(initial_cond)\n",
    "    cnt = 0; time = t0\n",
    "    sample_int = int(floor(tend / dt / nsamples)+1)\n",
    "    gfut = GridFunction(gfu.space,multidim=0)\n",
    "    gfut.AddMultiDimComponent(gfu.vec)\n",
    "    while time < tend - 0.5 * dt:\n",
    "        res = dt * f.vec - dt * a.mat * gfu.vec\n",
    "        gfu.vec.data += invmstar * res\n",
    "        print(\"\\r\",time,end=\"\")\n",
    "        scene.Redraw()\n",
    "        if cnt % sample_int == 0:\n",
    "            gfut.AddMultiDimComponent(gfu.vec)\n",
    "        cnt += 1; time = cnt * dt\n",
    "    return gfut"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "gfut = TimeStepping(invmstar, (1-y*y)*x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw(gfut, mesh, interpolate_multidim=True, animate=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Supplementary material:\n",
    "\n",
    "The simple time stepping can be extended in several directions:\n",
    "\n",
    "1. Instead of direct one could use iterative solvers\n",
    "2. R.h.s. data time-dependent\n",
    "3. boundary data time-dependent\n",
    "4. Replace implicit Euler time stepping with (Singly diagonal) Runge-Kutta methods\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Supplementary 1: Using iterative solvers\n",
    "\n",
    "* For a factorization of $M^\\ast$ (and ${M^\\ast}^{-1}$) we required a sparse matrix $M^\\ast$ \n",
    "* To store $M^\\ast$ **as a sparse matrix** requires new storage (and same memory layout of $A$ and $M$)\n",
    "* For iterative solvers we only require the matrix (and preconditioner) applications and it **suffices to have $M^\\ast$ available as a linear operator**\n",
    "* `mstar_alt = m.mat + dt * a.mat` has no storage but defines the operator action: sum of two matrix-vector multiplications\n",
    "\n",
    "iterative solver version (with Jacobi preconditining):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "mstar_alt = m.mat + dt * a.mat\n",
    "premstar_alt = m.mat.CreateSmoother() + dt * a.mat.CreateSmoother()\n",
    "print(premstar_alt)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now, replace the action of the inverse Matrix from the previous script with a solution of a conjugate gradient method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "from ngsolve.krylovspace import CGSolver\n",
    "invmstar_alt = CGSolver(mstar_alt, premstar_alt, tol=1e-8, printrates=\"\\r\", maxiter=200)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now, we can repeat the time stepping. As in 2D direct solvers are quite efficient, this simple preconditioned solver can hardly compete. We hence only do a few time steps:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "gfut_a1 = TimeStepping(invmstar_alt, (1-y*y)*x, tend=25*dt, nsamples=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "Draw(gfut_a1, mesh, interpolate_multidim=True, animate=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Supplementary 2: time-dependent r.h.s. data <a class=\"anchor\" id=\"TDRHS\"></a>\n",
    "\n",
    "Next: time-dependent r.h.s. data $f=f(t)$:\n",
    "\n",
    "* Use `Parameter` t representing the time. \n",
    "* A `Parameter` is a constant `CoefficientFunction` the value of which can be changed with the `Set`-function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "t = Parameter(0.0)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "An example of a time-dependent coefficient that we want to use as r.h.s. in the following is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "omega=1\n",
    "gausspt = exp(-6*((x+sin(omega*t))*(x+sin(omega*t))+y*y))-exp(-6*((x-sin(omega*t))*(x-sin(omega*t))+y*y))\n",
    "gff = GridFunction(L2(mesh,order=gfu.space.globalorder+1))\n",
    "gfft = GridFunction(gff.space,multidim=0)\n",
    "time = 0.0\n",
    "for i in range(7):\n",
    "    t.Set(3*i/6)\n",
    "    gff.Set(gausspt)\n",
    "    gfft.AddMultiDimComponent(gff.vec)\n",
    "Draw(gfft, mesh, interpolate_multidim=True, animate=True, autoscale=False, min=-1, max=1)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Accordingly, we define a different linear form which then has to be assembled in every time step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "ft = LinearForm(fes)\n",
    "ft += gausspt*v*dx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def TimeStepping_app2(invmstar, initial_cond = None, t0 = 0, tend = 2, \n",
    "                      nsamples = 10):\n",
    "    if initial_cond:\n",
    "        gfu.Set(initial_cond)\n",
    "    cnt = 0; time = t0\n",
    "    sample_int = int(floor(tend / dt / nsamples)+1)\n",
    "    gfut = GridFunction(gfu.space,multidim=0)\n",
    "    gfut.AddMultiDimComponent(gfu.vec)\n",
    "    while time < tend - 0.5 * dt:\n",
    "        t.Set(time)\n",
    "        ft.Assemble()\n",
    "        res = dt * ft.vec - dt * a.mat * gfu.vec\n",
    "        gfu.vec.data += invmstar * res\n",
    "        print(\"\\r\",time,end=\"\")\n",
    "        if cnt % sample_int == 0:\n",
    "            gfut.AddMultiDimComponent(gfu.vec)\n",
    "        cnt += 1; time = cnt * dt\n",
    "    return gfut"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "gfut_a2 = TimeStepping_app2(invmstar, initial_cond=CF(0),tend=2)\n",
    "Draw(gfut_a2, mesh, interpolate_multidim=True, animate=True, \n",
    "     min=-0.75,max=0.75,autoscale=False)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Supplementary 3: Time dependent boundary conditions <a class=\"anchor\" id=\"TDBND\"></a>\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "* $u|_{\\partial \\Omega} = u_D(t)$, $f=0$\n",
    "* implicit Euler time stepping method, non-incremental form:\n",
    "\n",
    "  $$\n",
    "    M^\\ast u^{n+1} = (M + \\Delta t A) u^{n+1} = M u^n\n",
    "  $$  \n",
    "  \n",
    "* Homogenize w.r.t. to boundary conditions, i.e. we split \n",
    "\n",
    "  $$ u^{n+1} = u^{n+1}_0 + u^{n+1}_D $$ \n",
    "  \n",
    "  where $u^{n+1}_D$ is a (discrete) function with correct boundary condition:  \n",
    "  \n",
    "$$\n",
    "  {M^\\ast} u^{n+1}_0 = M u^n - {M^\\ast} u^{n+1}_D\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "uD = CoefficientFunction( [(1-x*x)*IfPos(sin(0.3*pi*t),sin(0.3*pi*t),0),0,0,0])\n",
    "time = 0.0\n",
    "t.Set(0.0)\n",
    "gfu.Set(uD,BND)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     2,
     15
    ],
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def TimeStepping_app3(invmstar, initial_cond = None, t0 = 0, tend = 2, \n",
    "                      nsamples = 10):\n",
    "    if initial_cond:\n",
    "        gfu.Set(initial_cond)\n",
    "    cnt = 0; time = t0\n",
    "    sample_int = int(floor(tend / dt / nsamples)+1)\n",
    "    gfuD = GridFunction(gfu.space)\n",
    "    gfut = GridFunction(gfu.space,multidim=0)\n",
    "    gfut.AddMultiDimComponent(gfu.vec)\n",
    "    while time < tend - 0.5 * dt:\n",
    "        t.Set(time)\n",
    "        gfuD.Set(uD,BND)\n",
    "        res = m.mat * gfu.vec - mstar *gfuD.vec\n",
    "        gfu.vec.data = gfuD.vec + invmstar * res\n",
    "        print(\"\\r\",time,end=\"\")\n",
    "        if cnt % sample_int == 0:\n",
    "            gfut.AddMultiDimComponent(gfu.vec)\n",
    "        cnt += 1; time = cnt * dt\n",
    "    return gfut"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "gfut_a3 = TimeStepping_app3(invmstar, initial_cond=CF(0),tend=10)\n",
    "Draw(gfut_a3, mesh, interpolate_multidim=True, animate=True, \n",
    "     #settings = {\"subdivision\" : 10}, \n",
    "     deformation = True, min = 0, max = 0.2, autoscale = False)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Supplementary 4: Singly diagonally implicit Runge-Kutta methods <a class=\"anchor\" id=\"RK\"></a>\n",
    "\n",
    "\n",
    "Let us consider more sophisticated time integration methods, SDIRK methods. To simplify presentation we set $f=0$.\n",
    "\n",
    "SDIRK methods for the semi-discrete problem $\\frac{d}{dt} u = M^{-1} F(u) = -M^{-1} \\cdot A u$ are of the form:\n",
    "\n",
    "$$\n",
    "  u^{n+1} = u^n + \\Delta t M^{-1} \\sum_{i=0}^{s-1} b_i k_i\n",
    "$$\n",
    "\n",
    "with\n",
    "\n",
    "$$\n",
    "k_i = - A u_i \\quad \\text{ where $u_i$ is the solution to } \\quad \n",
    "(M + a_{ii} \\Delta t A) u_i = M u^n + \\Delta t \\sum_{j=0}^{i-1} a_{ij} k_j, \\quad i=0,..,s-1.\n",
    "$$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The coefficients a,b and c are stored in the butcher tableau:\n",
    "\n",
    "$$\n",
    "\\begin{array}{c|cccc}\n",
    "c_0 & a_{00} & 0 & \\ddots& \\\\\n",
    "c_1 & a_{10} & a_{11} & 0 & \\ddots \\\\\n",
    "\\vdots & \\vdots & \\vdots & \\ddots & 0\\\\\n",
    "c_{s-1} & a_{s-1,0} & a_{s-1,1} & \\dots& a_{s-1,s-1} \\\\ \\hline\n",
    " & b_{0} & b_1 & \\dots& b_{s-1} \\\\\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "For an SDIRK method we have $a^\\ast = a_{ii},~i=0,..,s-1$."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Simplest example: Implicit Euler\n",
    "$$\n",
    "\\begin{array}{c|c}\n",
    "1 & 1 \\\\ \\hline\n",
    " & 1  \\\\\n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "class sdirk1: #order 1 (implicit Euler)\n",
    "    stages = 1\n",
    "    a = [[1]]\n",
    "    b = [1]\n",
    "    c = [1]\n",
    "    astar = 1"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Midpoint rule:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "class MR: #order 2 (midpoint rule)\n",
    "    stages = 1\n",
    "    a = [[0.5]]\n",
    "    b = [1]\n",
    "    c = [0.5]\n",
    "    astar = 0.5"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can use for example the 2 stage SDIRK (order 3) method\n",
    "\n",
    "$$\n",
    "\\begin{array}{c|cc}\n",
    "p & p & 0 \\\\\n",
    "1-p & 1-2p & p \\\\ \\hline\n",
    " & 1/2 & 1/2 \\\\\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "with $p = \\frac{3 - \\sqrt{3}}{6}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "class sdirk2: #order 3\n",
    "    p = (3 - sqrt(3))/6\n",
    "    stages = 2\n",
    "    a = [[p, 0], \n",
    "         [1 - 2*p, p]]\n",
    "    b = [1/2, 1/2]\n",
    "    c = [p, 1 - p]\n",
    "    astar = p\n",
    "    "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "or the 5 stage SDIRK (order 4) method\n",
    "$$\n",
    "\\begin{array}{c|cccc}\n",
    "1/4 & 1/4  \\\\\n",
    "3/4 & 1/2 & 1/4 & \\\\\n",
    "11/20 & 17/50 & -1/25 & 1/4 \\\\\n",
    "1/2 & 371/1360 & -137/2720 & 15/544 & 1/4 \\\\\n",
    "1 & 25/4 & -49/48 & 125/16 & -85/12 & 1/4 \\\\ \\hline\n",
    " & 25/4 & -49/48 & 125/16 & -85/12 & 1/4 \\\\\n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "class sdirk5: #order 4\n",
    "    stages = 5\n",
    "    a = [[1/4, 0, 0, 0, 0], \n",
    "         [1/2, 1/4, 0, 0, 0], \n",
    "         [17/50,-1/25, 1/4, 0, 0],\n",
    "         [371/1360, -137/2720, 15/544, 1/4,0],\n",
    "         [25/24, -49/48, 125/16, -85/12, 1/4]]\n",
    "    b = [25/24, -49/48, 125/16, -85/12, 1/4]\n",
    "    c = [1/4, 3/4, 11/20, 1/2, 1]\n",
    "    astar = 1/4    "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We have to update the $M^\\ast$ matrix due to the different diagonal in the butcher tableau.\n",
    "Further, we will also need the inverse mass matrix to compute the slope states $k_i$. \n",
    "\n",
    "In the next function, we put setup of the matrices, their inverses and the time stepping in one (larger) loop:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def TimeStepping_app4(m, a, butchertab, dt,\n",
    "                      initial_cond = None, t0 = 0, tend = 2, \n",
    "                      nsamples = 10):\n",
    "    mstar.AsVector().data = m.mat.AsVector() + butchertab.astar * dt * a.mat.AsVector()\n",
    "    invmstar = mstar.Inverse(freedofs=fes.FreeDofs())\n",
    "    invmass = m.mat.Inverse(freedofs=fes.FreeDofs())    \n",
    "    \n",
    "    rhsi, Mu0, ui = gfu.vec.CreateVector(),gfu.vec.CreateVector(),gfu.vec.CreateVector()   \n",
    "    k = [gfu.vec.CreateVector() for i in range(butchertab.stages)]    \n",
    "\n",
    "    if initial_cond:\n",
    "        gfu.Set(initial_cond)\n",
    "    cnt = 0; time = t0\n",
    "    sample_int = int(floor(tend / dt / nsamples)+1)\n",
    "    gfuD = GridFunction(gfu.space)\n",
    "    gfut = GridFunction(gfu.space,multidim=0)\n",
    "    gfut.AddMultiDimComponent(gfu.vec)\n",
    "    while time < tend - 0.5 * dt:\n",
    "        Mu0.data = m.mat * gfu.vec\n",
    "        for i in range(butchertab.stages):\n",
    "            rhsi.data = Mu0\n",
    "            for j in range(0,i):\n",
    "                rhsi.data += dt * butchertab.a[i][j] * k[j]\n",
    "            t.Set(time+butchertab.c[i]*dt)\n",
    "            gfu.Set(uD,BND)\n",
    "            ui.data = gfu.vec; rhsi.data -= mstar * ui\n",
    "            ui.data += invmstar * rhsi\n",
    "            k[i].data = - a.mat * ui\n",
    "        t.Set(time+dt)\n",
    "        gfu.Set(uD,BND)\n",
    "        Mu0.data -= m.mat * gfu.vec\n",
    "        for i in range(0,butchertab.stages):\n",
    "            Mu0.data += dt * butchertab.b[i] * k[i]\n",
    "        gfu.vec.data += invmass * Mu0        \n",
    "        print(\"\\r\",time,end=\"\")\n",
    "        if cnt % sample_int == 0:\n",
    "            gfut.AddMultiDimComponent(gfu.vec)\n",
    "        cnt += 1; time = cnt * dt\n",
    "    return gfut"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Next, we simply apply the time stepping:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "gfut_a4 = TimeStepping_app4(m, a, butchertab = sdirk5(), \n",
    "                            dt = 0.25, initial_cond=CF(0),tend=10)\n",
    "Draw(gfut_a4, mesh, interpolate_multidim=True, animate=True, \n",
    "     order = 3, \n",
    "     deformation = True, min = 0, max = 0.2, autoscale = False)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Tasks**\n",
    "\n",
    "If you want to play around and explore this tutorial unit, here are a few suggestions:\n",
    "\n",
    "* Make a convergence study (for convergence in space, time and space and time)\n",
    "* Re-introduce non-zero r.h.s. in the RK-script\n",
    "* Generalize the RK script to fully implicit RK schemes (using product spaces, or block matrices, ...)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "interpreter": {
   "hash": "767d51c1340bd893661ea55ea3124f6de3c7a262a8b4abca0554b478b1e2ff90"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.10"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}