{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# 3.7 DG/HDG splitting\n",
    "\n",
    "When solving unsteady problems with an operator-splitting method it might be beneficial to consider different space discretizations for different operators.\n",
    "\n",
    "For a problem of the form\n",
    "\n",
    "$$\n",
    "  \\partial_t u + A u + C u = 0\n",
    "$$\n",
    "\n",
    "We consider the operator splitting:\n",
    "\n",
    "* 1.Step: Given $u^0$, solve $t^n \\to t^{n+1}$: $\\quad \\partial_t u + C u = 0 \\Rightarrow u^{\\frac12}$\n",
    "* 2.Step: Given $u^{\\frac12}$, solve $t^n \\to t^{n+1}$: $\\quad \\partial_t u + A u = 0 \\Rightarrow u^{1}$\n",
    "\n",
    "This splitting is only first order accurate but allows to take different time discretizations for the substeps 1 and 2. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In this example we consider the Navier-Stokes problem again (cf. [unit 3.2](../unit-3.2-navierstokes/navierstokes.ipynb)) and want to combine \n",
    "\n",
    "* an $H(div)$ -conforming Hybrid DG method (which is a very good discretization for Stokes-type problems) and\n",
    "* a standard **upwind DG** method for the discretization of the convection\n",
    "\n",
    "The weak form is: Find $(\\mathbf{u},p):[0,T] \\to (H_{0,D}^1)^d \\times L^2$, s.t.\n",
    "\n",
    "\\begin{align}\n",
    "\\int_{\\Omega} \\partial_t \\mathbf{u} \\cdot v + \\int_{\\Omega} \\nu \\nabla \\mathbf{u} \\nabla \\mathbf{v} + \\mathbf{u} \\cdot \\nabla \\mathbf{u} \\ \\mathbf{v} - \\int_{\\Omega} \\operatorname{div}(\\mathbf{v}) p &= \\int f \\mathbf{v}  && \\forall \\mathbf{v} \\in (H_{0,D}^1)^d, \\\\ \n",
    "- \\int_{\\Omega} \\operatorname{div}(\\mathbf{u}) q &= 0 && \\forall q \\in L^2, \\\\\n",
    "\\quad \\mathbf{u}(t=0) & = \\mathbf{u}_0\n",
    "\\end{align}\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Again, we consider the benchmark setup from http://www.featflow.de/en/benchmarks/cfdbenchmarking/flow/dfg_benchmark2_re100.html . The geometry:\n",
    "\n",
    "![](../unit-3.2-navierstokes/geometry.png)\n",
    "\n",
    "The viscosity is set to $\\nu = 10^{-3}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "order = 3\n",
    "# import libraries, set up geometry and generate mesh\n",
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "from netgen.occ import *\n",
    "shape = Rectangle(2,0.41).Circle(0.2,0.2,0.05).Reverse().Face()\n",
    "shape.edges.name=\"cyl\"\n",
    "shape.edges.Min(X).name=\"inlet\"\n",
    "shape.edges.Max(X).name=\"outlet\"\n",
    "shape.edges.Min(Y).name=\"wall\"\n",
    "shape.edges.Max(Y).name=\"wall\"\n",
    "mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.07)).Curve(order)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For the HDG formulation we use the product space with\n",
    "\n",
    "* $BDM_k$: $H(div)$ conforming FE space (degree k)\n",
    "* Vector facet space: facet functions of degree k (vector valued and only in tangential direction)\n",
    "* piecewise polynomials up to degree $k-1$ for the pressure\n",
    "\n",
    "![](resources/stokeshdghdiv.png)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## HDG spaces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "V1 = HDiv ( mesh, order = order, dirichlet = \"wall|cyl|inlet\" )\n",
    "V2 = TangentialFacetFESpace(mesh, order = order, dirichlet = \"wall|cyl|inlet\" )\n",
    "Q = L2( mesh, order = order-1)\n",
    "V = FESpace ([V1,V2,Q])\n",
    "\n",
    "u, uhat, p = V.TrialFunction()\n",
    "v, vhat, q = V.TestFunction()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parameter and directions (normal / tangential projection):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "nu = 0.001 # viscosity\n",
    "n = specialcf.normal(mesh.dim)\n",
    "h = specialcf.mesh_size\n",
    "def tang(vec):\n",
    "    return vec - (vec*n)*n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Stokes discretization\n",
    "The bilinear form to the HDG (see [unit-2.8-DG](../unit-2.8-DG/DG.ipynb) for scalar HDG) discretized Stokes problem is:\n",
    "\n",
    "\\begin{align*}\n",
    "K_h((\\mathbf{u}_T,\\mathbf{u}_F,p),(\\mathbf{v}_T,\\mathbf{v}_F,q) := & \\displaystyle \\sum_{T\\in\\mathcal{T}} \\left\\{ \\int_{T} \\nu {\\nabla} {\\mathbf{u}_T} \\! : \\! {\\nabla} {\\mathbf{v}_T} \\ d {x} - \\int_{\\Omega} \\operatorname{div}(\\mathbf{u}) q \\ d {x} - \\int_{\\Omega} \\operatorname{div}(\\mathbf{v}) p \\ d {x} \\right. \\\\\n",
    "& \\qquad \\left. - \\int_{\\partial T} \\nu \\frac{\\partial {\\mathbf{u}_T}}{\\partial {n} }  [ \\mathbf{v} ]_t \\ d {s}\n",
    "- \\int_{\\partial T} \\nu \\frac{\\partial {\\mathbf{v}_T}}{\\partial {n} }  [ \\mathbf{u} ]_t \\ d {s}\n",
    "+ \\int_{\\partial T} \\nu \\frac{\\alpha k^2}{h} [\\mathbf{u}]_t [\\mathbf{v}]_t \\ d {s}  \\right\\}\n",
    "\\end{align*}\n",
    "\n",
    "where $[\\cdot]_t$ is the tangential projection of the jump $(\\cdot)_T - (\\cdot)_F$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "alpha = 4  # SIP parameter\n",
    "dS = dx(element_boundary=True)\n",
    "a = BilinearForm ( V, symmetric=True)\n",
    "a += nu * InnerProduct ( Grad(u), Grad(v) ) * dx\n",
    "a += nu * InnerProduct ( Grad(u)*n, tang(vhat-v) ) * dS \n",
    "a += nu * InnerProduct ( Grad(v)*n, tang(uhat-u) ) * dS\n",
    "a += nu * alpha*order*order/h * InnerProduct(tang(vhat-v), tang(uhat-u)) * dS\n",
    "a += (-div(u)*q - div(v)*p) * dx\n",
    "a.Assemble()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The mass matrix is simply\n",
    "\n",
    "\\begin{align*}\n",
    "M_h((\\mathbf{u}_T,\\mathbf{u}_F,p),(\\mathbf{v}_T,\\mathbf{v}_F,q) := \\int_{\\Omega} \\mathbf{u}_T \\cdot \\mathbf{v}_T dx\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "m = BilinearForm(u * v * dx, symmetric=True).Assemble()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Initial conditions\n",
    "As initial condition we solve the Stokes problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "gfu = GridFunction(V)\n",
    "\n",
    "U0 = 1.5\n",
    "uin = CF( (U0*4*y*(0.41-y)/(0.41*0.41),0) )\n",
    "gfu.components[0].Set(uin, definedon=mesh.Boundaries(\"inlet\"))\n",
    "\n",
    "invstokes = a.mat.Inverse(V.FreeDofs(), inverse=\"umfpack\")\n",
    "gfu.vec.data += invstokes @ -a.mat * gfu.vec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "scenes = [ \\\n",
    "    Draw( gfu.components[0], mesh, \"velocity\"), \\\n",
    "    Draw( gfu.components[2], mesh, \"pressure\")];"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Application of convection\n",
    "\n",
    "In the operator splitted approach we want to apply only operator applications for the convection part.\n",
    "Further, we want to do this in a usual DG setting. \n",
    "As a model problem we use the following procedure:\n",
    "\n",
    "* Given $(\\mathbf{u},p)$ in HDG space: project into $\\hat{\\mathbf{u}}$ in usual DG space\n",
    "* Solve $\\partial_t \\hat{\\mathbf{u}} = C \\hat{\\mathbf{u}}$ by explicit scheme (involves convection evaluations and mass matrix operations only)\n",
    "* Solve Unsteady Stokes step with r.h.s. from convection sub-problem. To this end evaluate mixed mass matrix $\\int \\hat{\\mathbf{u}} \\cdot \\mathbf{v}$ to obtain a functional on the HDG space"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For the projection steps we use mixed mass matrices:\n",
    "\n",
    "* $M_m$ : $HDG \\times DG \\to \\mathbb{R}$\n",
    "* $M_m^T$ : $DG \\times HDG \\to \\mathbb{R}$\n",
    "* $M_{DG}$ : $DG \\times DG \\to \\mathbb{R}$ (block diagonal)\n",
    "\n",
    "To set up mixed mass matrices we use a bilinear form with two different FESpaces. \n",
    "\n",
    "We do not assemble the matrices as we will only need the matrix-vector applications of $M_m$, $M_m^T$ and $M_{DG}^{-1}$.\n",
    "\n",
    "There is a more elegant version (`ConvertL2Operator`, similar to `TraceOperator`) that we discuss [below](#convertl2)."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### mixed mass matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "VL2 = VectorL2(mesh, order=order, piola=True)\n",
    "uL2, vL2 = VL2.TnT()\n",
    "bfmixed = BilinearForm ( vL2*u*dx, nonassemble=True )\n",
    "bfmixedT = BilinearForm ( uL2*v*dx, nonassemble=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### convection operator\n",
    "* convection operation with standard Upwinding\n",
    "* No set up of the matrix (only interested in operator applications)\n",
    "* for the advection velocity we use the $H(div)$-conforming velocity (which is pointwise divergence free)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "vel = gfu.components[0]\n",
    "convL2 = BilinearForm( (-InnerProduct(Grad(vL2) * vel, uL2)) * dx, nonassemble=True )\n",
    "un = InnerProduct(vel,n)\n",
    "upwindL2 = IfPos(un, un*uL2, un*uL2.Other(bnd=uin))\n",
    "\n",
    "dskel_inner  = dx(skeleton=True)\n",
    "dskel_bound  = ds(skeleton=True)\n",
    "\n",
    "convL2 += InnerProduct (upwindL2, vL2-vL2.Other()) * dskel_inner\n",
    "convL2 += InnerProduct (upwindL2, vL2) * dskel_bound"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Solution of convection steps as an operator (`BaseMatrix`)\n",
    "We now define the solution of the convection problem for an initial data $u$ in the HDG space. The return value (`y`) is $M_m \\hat{u}$ where $\\hat{u}$ is the solution of several explicit Euler steps of the convection problem\n",
    "$$\n",
    "  \\partial_t \\hat{u} = - M_{DG}^{-1} C \\hat{u}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "class PropagateConvection(BaseMatrix):\n",
    "    def __init__(self,tau,steps):\n",
    "        super(PropagateConvection, self).__init__()\n",
    "        self.tau = tau; self.steps = steps\n",
    "        self.h = V.ndof; self.w = V.ndof # operator domain and range\n",
    "        self.mL2 = VL2.Mass(Id(mesh.dim)); self.invmL2 = self.mL2.Inverse()\n",
    "        self.vecL2 = bfmixed.mat.CreateColVector() # temp vector\n",
    "    def Mult(self, x, y):\n",
    "        self.vecL2.data = self.invmL2 @ bfmixed.mat * x # <- project from Hdiv to L2\n",
    "        for i in range(self.steps):\n",
    "            self.vecL2.data -= self.tau/self.steps * self.invmL2 @ convL2.mat * self.vecL2\n",
    "        y.data = bfmixedT.mat * self.vecL2\n",
    "    def CreateColVector(self):\n",
    "        return CreateVVector(self.h)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The operator splitted method (Yanenko splitting)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Setup of $M^*$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "t = 0; tend = 0\n",
    "tau = 0.01; substeps = 10\n",
    "\n",
    "mstar = m.mat.CreateMatrix()\n",
    "mstar.AsVector().data = m.mat.AsVector() + tau * a.mat.AsVector()\n",
    "inv = mstar.Inverse(V.FreeDofs(), inverse=\"umfpack\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The time loop:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "for s in scenes : s.Draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "tend += 1\n",
    "res = gfu.vec.CreateVector()\n",
    "convpropagator = PropagateConvection(tau,substeps)\n",
    "while t < tend:\n",
    "    gfu.vec.data += inv @ (convpropagator - mstar) * gfu.vec\n",
    "    t += tau\n",
    "    print (\"\\r  t =\", t, end=\"\")\n",
    "    for s in scenes : s.Redraw()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## `ConvertL2Operator` <a class=\"anchor\" id=\"convertl2\"></a>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Instead of doing the conversion between the spaces \"by hand\", we can use the convenience operator `ConvertL2Operator` that realizes\n",
    "$ M_{DG}^{-1} \\cdot M_m $:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "HDG_to_L2 = V1.ConvertL2Operator(VL2) @ Embedding(V.ndof, V.Range(0)).T"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "The transpose of $M_{DG}^{-1} \\cdot M_m$ is $M_m^T \\cdot M_{DG}^{-1}$\n",
    "* `V1.ConvertL2Operator(VL2)`: V1 $\\to$ VL2\n",
    "* `V1.ConvertL2Operator(VL2).T`: VL2' $\\to$ V1'"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Overwrite the application and use the `ConvertL2Operator` instead ot the mixed mass matrices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "def NewMult(self, x, y):\n",
    "    self.vecL2.data = HDG_to_L2 * x # <- project from Hdiv to L2\n",
    "    for i in range(self.steps):\n",
    "        self.vecL2.data -= tau/self.steps*self.invmL2@convL2.mat*self.vecL2\n",
    "    y.data = HDG_to_L2.T @ self.mL2 * self.vecL2    \n",
    "    \n",
    "PropagateConvection.Mult = NewMult"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The time loop:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "for s in scenes : s.Draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "tend += 1\n",
    "convpropagator = PropagateConvection(tau,substeps)\n",
    "while t < tend:\n",
    "    gfu.vec.data += inv @ (convpropagator - mstar) * gfu.vec\n",
    "    t += tau\n",
    "    print (\"\\r  t =\", t, end=\"\")\n",
    "    for s in scenes : s.Redraw()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Tasks**:\n",
    "\n",
    "* Implement higher order Splitting (e.g. Strang) schemes\n",
    "* Implement IMEX schemes and compare"
   ]
  }
 ],
 "metadata": {
  "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
}