NGSolve24

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# 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": 57,
   "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": 58,
   "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": 59,
   "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": 60,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3dab9f6689fe4457bb2a19981a74a835",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.2…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "BaseWebGuiScene"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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": 61,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<ngsolve.comp.BilinearForm at 0x1f3ac9b7cf0>"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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": 62,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "m.mat.nze = 11257, a.mat.nze=11257, mstar.nze=11257\n"
     ]
    }
   ],
   "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": 63,
   "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": 64,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4b5a72709cf04042b1da99f095d02a88",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.2…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<ngsolve.comp.LinearForm at 0x1f380aad870>"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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": 65,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "713c47beb752412489700873f09514bb",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.2…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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": 66,
   "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": 67,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 0.143000000000000024"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "File \u001b[1;32m<timed exec>:1\u001b[0m\n",
      "Cell \u001b[1;32mIn[66], line 15\u001b[0m, in \u001b[0;36mTimeStepping\u001b[1;34m(invmstar, initial_cond, t0, tend, nsamples)\u001b[0m\n\u001b[0;32m     13\u001b[0m gfu\u001b[38;5;241m.\u001b[39mvec\u001b[38;5;241m.\u001b[39mdata \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m invmstar \u001b[38;5;241m*\u001b[39m res\n\u001b[0;32m     14\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;130;01m\\r\u001b[39;00m\u001b[38;5;124m\"\u001b[39m,time,end\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m---> 15\u001b[0m \u001b[43mscene\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mRedraw\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m     16\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m cnt \u001b[38;5;241m%\u001b[39m sample_int \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m0\u001b[39m:\n\u001b[0;32m     17\u001b[0m     gfut\u001b[38;5;241m.\u001b[39mAddMultiDimComponent(gfu\u001b[38;5;241m.\u001b[39mvec)\n",
      "File \u001b[1;32mc:\\Users\\User\\OneDrive\\Desktop\\ngs24\\Lib\\site-packages\\netgen\\webgui.py:234\u001b[0m, in \u001b[0;36mWebGLScene.Redraw\u001b[1;34m(self, *args, **kwargs)\u001b[0m\n\u001b[0;32m    232\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39margs \u001b[38;5;241m=\u001b[39m new_scene\u001b[38;5;241m.\u001b[39margs\n\u001b[0;32m    233\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mkwargs \u001b[38;5;241m=\u001b[39m new_scene\u001b[38;5;241m.\u001b[39mkwargs\n\u001b[1;32m--> 234\u001b[0m \u001b[38;5;28;43msuper\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mRedraw\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[1;32mc:\\Users\\User\\OneDrive\\Desktop\\ngs24\\Lib\\site-packages\\webgui_jupyter_widgets\\widget.py:81\u001b[0m, in \u001b[0;36mBaseWebGuiScene.Redraw\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m     79\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mRedraw\u001b[39m(\u001b[38;5;28mself\u001b[39m):\n\u001b[0;32m     80\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mencoding\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mbinary\u001b[39m\u001b[38;5;124m'\u001b[39m\n\u001b[1;32m---> 81\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mwidget\u001b[38;5;241m.\u001b[39mvalue \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mGetData\u001b[49m\u001b[43m(\u001b[49m\u001b[43mset_minmax\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mFalse\u001b[39;49;00m\u001b[43m)\u001b[49m\n",
      "File \u001b[1;32mc:\\Users\\User\\OneDrive\\Desktop\\ngs24\\Lib\\site-packages\\netgen\\webgui.py:241\u001b[0m, in \u001b[0;36mWebGLScene.GetData\u001b[1;34m(self, set_minmax)\u001b[0m\n\u001b[0;32m    239\u001b[0m d \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mNone\u001b[39;00m\n\u001b[0;32m    240\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mtype\u001b[39m(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mobj) \u001b[38;5;129;01min\u001b[39;00m _registered_draw_types:\n\u001b[1;32m--> 241\u001b[0m     d \u001b[38;5;241m=\u001b[39m \u001b[43m_registered_draw_types\u001b[49m\u001b[43m[\u001b[49m\u001b[43mtyp\u001b[49m\u001b[43m]\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mobj\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m    242\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m    243\u001b[0m     \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01minspect\u001b[39;00m\n",
      "File \u001b[1;32mc:\\Users\\User\\OneDrive\\Desktop\\ngs24\\Lib\\site-packages\\ngsolve\\webgui.py:162\u001b[0m, in \u001b[0;36mGetData\u001b[1;34m(obj, args, kwargs)\u001b[0m\n\u001b[0;32m    160\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mintpoints\u001b[39m\u001b[38;5;124m'\u001b[39m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;129;01min\u001b[39;00m kwargs:\n\u001b[0;32m    161\u001b[0m     kwargs[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mintpoints\u001b[39m\u001b[38;5;124m'\u001b[39m]  \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mNone\u001b[39;00m\n\u001b[1;32m--> 162\u001b[0m d \u001b[38;5;241m=\u001b[39m \u001b[43mBuildRenderData\u001b[49m\u001b[43m(\u001b[49m\u001b[43mmesh\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcf\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43morder\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43morder\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdraw_surf\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mdraw_surf\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdraw_vol\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mdraw_vol\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mintpoints\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mintpoints\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdeformation\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mdeformation\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mregions\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mregions\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mobjects\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mobjects\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mnodal_p1\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mnodal_p1\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43msettings\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43msettings\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m    164\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(cf, ngs\u001b[38;5;241m.\u001b[39mGridFunction) \u001b[38;5;129;01mand\u001b[39;00m \u001b[38;5;28mlen\u001b[39m(cf\u001b[38;5;241m.\u001b[39mvecs)\u001b[38;5;241m>\u001b[39m\u001b[38;5;241m1\u001b[39m:\n\u001b[0;32m    165\u001b[0m     \u001b[38;5;66;03m# multidim gridfunction - generate data for each component\u001b[39;00m\n\u001b[0;32m    166\u001b[0m     gf \u001b[38;5;241m=\u001b[39m ngs\u001b[38;5;241m.\u001b[39mGridFunction(cf\u001b[38;5;241m.\u001b[39mspace)\n",
      "File \u001b[1;32mc:\\Users\\User\\OneDrive\\Desktop\\ngs24\\Lib\\site-packages\\ngsolve\\webgui.py:406\u001b[0m, in \u001b[0;36mBuildRenderData\u001b[1;34m(mesh, func, order, draw_surf, draw_vol, intpoints, deformation, regions, objects, nodal_p1, encoding, settings)\u001b[0m\n\u001b[0;32m    404\u001b[0m timer2list\u001b[38;5;241m.\u001b[39mStart()        \n\u001b[0;32m    405\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m i \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(og\u001b[38;5;241m+\u001b[39m\u001b[38;5;241m1\u001b[39m):\n\u001b[1;32m--> 406\u001b[0m     Bezier_points\u001b[38;5;241m.\u001b[39mappend(\u001b[43mencode\u001b[49m\u001b[43m(\u001b[49m\u001b[43mBezierPnts\u001b[49m\u001b[43m[\u001b[49m\u001b[43mi\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdtype\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfloat32\u001b[49m\u001b[43m)\u001b[49m)\n\u001b[0;32m    407\u001b[0m timer2list\u001b[38;5;241m.\u001b[39mStop()        \n\u001b[0;32m    409\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m func2 \u001b[38;5;129;01mand\u001b[39;00m draw_surf:\n",
      "File \u001b[1;32mc:\\Users\\User\\OneDrive\\Desktop\\ngs24\\Lib\\site-packages\\ngsolve\\webgui.py:265\u001b[0m, in \u001b[0;36mBuildRenderData.<locals>.encode\u001b[1;34m(*args, **kwargs)\u001b[0m\n\u001b[0;32m    264\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mencode\u001b[39m(\u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs):\n\u001b[1;32m--> 265\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mencodeData\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mencoding\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mencoding\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[1;32mc:\\Users\\User\\OneDrive\\Desktop\\ngs24\\Lib\\site-packages\\netgen\\webgui.py:17\u001b[0m, in \u001b[0;36mencodeData\u001b[1;34m(data, dtype, encoding)\u001b[0m\n\u001b[0;32m     15\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mbase64\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m b64encode\n\u001b[0;32m     16\u001b[0m dtype \u001b[38;5;241m=\u001b[39m dtype \u001b[38;5;129;01mor\u001b[39;00m data\u001b[38;5;241m.\u001b[39mdtype\n\u001b[1;32m---> 17\u001b[0m values \u001b[38;5;241m=\u001b[39m \u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43marray\u001b[49m\u001b[43m(\u001b[49m\u001b[43mdata\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mflatten\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdtype\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdtype\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m     18\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m encoding\u001b[38;5;241m==\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mb64\u001b[39m\u001b[38;5;124m'\u001b[39m:\n\u001b[0;32m     19\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m b64encode(values)\u001b[38;5;241m.\u001b[39mdecode(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mascii\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "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": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cd1f0d1b3df04480bc2b331962d0e951",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.2…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "BaseWebGuiScene"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sum of\n",
      "Scale a = 1\n",
      "Print base-matrix\n",
      "Scale b = 1\n",
      "Scale with 0.001:\n",
      "Print base-matrix\n",
      "\n"
     ]
    }
   ],
   "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": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 0.989000000000000144"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "File \u001b[1;32m<timed exec>:1\u001b[0m\n",
      "Cell \u001b[1;32mIn[14], line 13\u001b[0m, in \u001b[0;36mTimeStepping\u001b[1;34m(invmstar, initial_cond, t0, tend, nsamples)\u001b[0m\n\u001b[0;32m     11\u001b[0m \u001b[38;5;28;01mwhile\u001b[39;00m time \u001b[38;5;241m<\u001b[39m tend \u001b[38;5;241m-\u001b[39m \u001b[38;5;241m0.5\u001b[39m \u001b[38;5;241m*\u001b[39m dt:\n\u001b[0;32m     12\u001b[0m     res \u001b[38;5;241m=\u001b[39m dt \u001b[38;5;241m*\u001b[39m f\u001b[38;5;241m.\u001b[39mvec \u001b[38;5;241m-\u001b[39m dt \u001b[38;5;241m*\u001b[39m a\u001b[38;5;241m.\u001b[39mmat \u001b[38;5;241m*\u001b[39m gfu\u001b[38;5;241m.\u001b[39mvec\n\u001b[1;32m---> 13\u001b[0m     \u001b[43mgfu\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mvec\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mdata\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m+\u001b[39;49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m \u001b[49m\u001b[43minvmstar\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mres\u001b[49m\n\u001b[0;32m     14\u001b[0m     \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;130;01m\\r\u001b[39;00m\u001b[38;5;124m\"\u001b[39m,time,end\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[0;32m     15\u001b[0m     scene\u001b[38;5;241m.\u001b[39mRedraw()\n",
      "File \u001b[1;32mc:\\Users\\User\\OneDrive\\Desktop\\ngs24\\Lib\\site-packages\\ngsolve\\krylovspace.py:113\u001b[0m, in \u001b[0;36mLinearSolver.Mult\u001b[1;34m(self, x, y)\u001b[0m\n\u001b[0;32m    112\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mMult\u001b[39m(\u001b[38;5;28mself\u001b[39m, x : BaseVector, y : BaseVector) \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m>\u001b[39m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m--> 113\u001b[0m     \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mSolve\u001b[49m\u001b[43m(\u001b[49m\u001b[43mrhs\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mx\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43msol\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43my\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43minitialize\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mTrue\u001b[39;49;00m\u001b[43m)\u001b[49m\n",
      "File \u001b[1;32mc:\\Users\\User\\OneDrive\\Desktop\\ngs24\\Lib\\site-packages\\netgen\\__init__.py:56\u001b[0m, in \u001b[0;36mTimeFunction.<locals>.retfunc\u001b[1;34m(*args, **kwargs)\u001b[0m\n\u001b[0;32m     54\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mretfunc\u001b[39m(\u001b[38;5;241m*\u001b[39margs,\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs):\n\u001b[0;32m     55\u001b[0m     \u001b[38;5;28;01mwith\u001b[39;00m timer:\n\u001b[1;32m---> 56\u001b[0m         ret \u001b[38;5;241m=\u001b[39m \u001b[43mfunc\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m     57\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m ret\n",
      "File \u001b[1;32mc:\\Users\\User\\OneDrive\\Desktop\\ngs24\\Lib\\site-packages\\ngsolve\\krylovspace.py:94\u001b[0m, in \u001b[0;36mLinearSolver.Solve\u001b[1;34m(self, rhs, sol, initialize)\u001b[0m\n\u001b[0;32m     92\u001b[0m     sol[:] \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m0\u001b[39m\n\u001b[0;32m     93\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39msol \u001b[38;5;241m=\u001b[39m sol\n\u001b[1;32m---> 94\u001b[0m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_SolveImpl\u001b[49m\u001b[43m(\u001b[49m\u001b[43mrhs\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mrhs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43msol\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43msol\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m     95\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m old_status[\u001b[38;5;241m0\u001b[39m] \u001b[38;5;241m!=\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124midle\u001b[39m\u001b[38;5;124m\"\u001b[39m:\n\u001b[0;32m     96\u001b[0m     _PushStatus(old_status[\u001b[38;5;241m0\u001b[39m])\n",
      "File \u001b[1;32mc:\\Users\\User\\OneDrive\\Desktop\\ngs24\\Lib\\site-packages\\ngsolve\\krylovspace.py:229\u001b[0m, in \u001b[0;36mCGSolver._SolveImpl\u001b[1;34m(self, rhs, sol)\u001b[0m\n\u001b[0;32m    227\u001b[0m w\u001b[38;5;241m.\u001b[39mdata \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mmat \u001b[38;5;241m*\u001b[39m s\n\u001b[0;32m    228\u001b[0m wd \u001b[38;5;241m=\u001b[39m wdn\n\u001b[1;32m--> 229\u001b[0m as_s \u001b[38;5;241m=\u001b[39m \u001b[43ms\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mInnerProduct\u001b[49m\u001b[43m(\u001b[49m\u001b[43mw\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mconjugate\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mconjugate\u001b[49m\u001b[43m)\u001b[49m        \n\u001b[0;32m    230\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m as_s \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m0\u001b[39m \u001b[38;5;129;01mor\u001b[39;00m wd \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m0\u001b[39m: \u001b[38;5;28;01mbreak\u001b[39;00m\n\u001b[0;32m    231\u001b[0m alpha \u001b[38;5;241m=\u001b[39m wd \u001b[38;5;241m/\u001b[39m as_s\n",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "%%time\n",
    "gfut_a1 = TimeStepping(invmstar_alt, (1-y*y)*x )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'gfut_a1' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[1;32mIn[20], line 1\u001b[0m\n\u001b[1;32m----> 1\u001b[0m Draw(\u001b[43mgfut_a1\u001b[49m, mesh, interpolate_multidim\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m, animate\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m)\n",
      "\u001b[1;31mNameError\u001b[0m: name 'gfut_a1' is not defined"
     ]
    }
   ],
   "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": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "56f0d0c1dc044d08bf78c1f4eb02884f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.2…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "BaseWebGuiScene"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 4.999000000000000544"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e735d5fb91b249c6a7a9ff7adb867d75",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.2…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: total: 3.92 s\n",
      "Wall time: 4.25 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "BaseWebGuiScene"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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": 85,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5acb90ee2f2a4a4b8823685d9f85a5fc",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.2…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "BaseWebGuiScene"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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": 86,
   "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": 87,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 1.999CPU times: total: 1.2 s\n",
      "Wall time: 1.48 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "gfut_a3 = TimeStepping_app3(invmstar, initial_cond=CF(gaussp),tend=2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7d1ee3c2bdbe4232a92f64aa27f20e08",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='500px', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.2…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "BaseWebGuiScene"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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
}