NGSolve24

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "(complementary:time)=\n",
    "# Time dependant problems\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}\")"
   ]
  },
  {
   "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",
    "\n",
    "# We cannot do this:\n",
    "# mstar.data = m.mat + dt * a.mat\n",
    "# unsupported operand type(s) for *: 'NoneType' and 'ngsolve.la.DynamicVectorExpression'\n",
    "\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",
    "    \n",
    "    gfut = GridFunction(gfu.space,multidim=0)\n",
    "    gfut.AddMultiDimComponent(gfu.vec)\n",
    "\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": "markdown",
   "metadata": {},
   "source": [
    " 1.999CPU times: user 7.05 s, sys: 630 ms, total: 7.68 s\n",
    "Wall time: 6.57 s"
   ]
  },
  {
   "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": [
    "## 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": "markdown",
   "metadata": {},
   "source": [
    "** does not perform well for this problem, just for explanation **\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "mstar_alt = m.mat + dt * a.mat\n",
    "\n",
    "precond = preconditioners.Local\n",
    "preflags ={\"GS\":False}\n",
    "\n",
    "prea = precond(a, **preflags)\n",
    "prem = precond(m, **preflags)\n",
    "\n",
    "\n",
    "premstar_alt = prem.mat + dt * prea.mat\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, plotrates=False, 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 )"
   ]
  },
  {
   "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": [
    "## time-dependent r.h.s. data \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=pi\n",
    "gausspt = exp(-10*((x-0.5)**2 + (y-0.5)**2) )*(sin(omega*t))\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, min=-1, max=1, autoscale=False);"
   ]
  },
  {
   "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 = 20):\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=5)\n",
    "Draw(gfut_a2, mesh, interpolate_multidim=True, animate=True, \n",
    "     min=-0.25,max=0.25,autoscale=False);"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Time dependent boundary conditions\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( (cos(5*(x+t))*y ))\n",
    "time = 0.0\n",
    "t.Set(0.0)\n",
    "gfu.Set(uD,BND)\n",
    "Draw(gfu,mesh,\"uD\");"
   ]
  },
  {
   "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(gaussp),tend=2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "Draw(gfut_a3, mesh, interpolate_multidim=True, animate=True, \n",
    "     #settings = {\"subdivision\" : 10}, \n",
    "     deformation = False, min = 0, max = 1, autoscale = True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.12.4"
  },
  "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
}