{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 8.4 Space-time discretizations on fitted geometry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial, we consider the heat equation with Dirichlet boundary conditions and want to solve it with a fitted DG-in-time space-time finite element method. The problem is posed on the domain $\\Omega = [0,1]^2$ and time in $[0,1]$. In detail, the PDE reads:\n",
    "\n",
    "$$\n",
    "\\left\\{\n",
    "\\begin{aligned}\n",
    "\\partial_t u - \\Delta u &= f \\quad \\text{ in } \\Omega \\text{ for all } t \\in [0,1],  & \\\\\n",
    "~ \\partial_{\\mathbf{n}} u &=  0  \\quad \\text{ on } \\partial \\Omega, & \\\\\n",
    "u &= u_0  \\quad \\text{at } t=0. & \\\\\n",
    "\\end{aligned}\\right.\n",
    "$$\n",
    "\n",
    "We calculate the right-hand side $f$ such that the solution is $u = \\sin(\\pi t) \\cdot \\sin(\\pi  x)^2 \\cdot \\sin(\\pi y)^2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.geom2d import unit_square\n",
    "from ngsolve import *\n",
    "from xfem import *\n",
    "import time\n",
    "from math import pi\n",
    "ngsglobals.msg_level = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We opt for a first-order method in space and time and choose an appropriate timestep."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Space finite element order\n",
    "order = 1\n",
    "# Time finite element order\n",
    "k_t = 1\n",
    "# Final simulation time\n",
    "tend = 1.0\n",
    "# Time step\n",
    "delta_t = 1 / 32"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In `ngsxfem` all space-time objects are designed for the reference time variable `tref` living on the reference time interval $[0,1]$. To obtain integrals and time derivatives with respect to the current time interval corresponding operators or objects need to be rescaled:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def dt(u):\n",
    "    return 1.0 / delta_t * dtref(u)\n",
    "\n",
    "tnew = 0\n",
    "told = Parameter(0)\n",
    "t = told + delta_t * tref    \n",
    "t.MakeVariable()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we opt for a standard space-time finite element discretisation based on the previous parameters. The space-time finite element space is of tensor-product form where the finite element for the time direction, a `ScalarTimeFE` uses a nodal basis. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = Mesh(unit_square.GenerateMesh(maxh=0.05, quad_dominated=False))\n",
    "\n",
    "V = H1(mesh, order=order, dirichlet=\".*\")\n",
    "tfe = ScalarTimeFE(k_t)\n",
    "st_fes = tfe * V # tensor product of a space finite element space and a 1D finite element"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For evaluation and time stepping the extraction of spatial functions from a space-time function will be required at several places. To obtain a spatial gridfunction from the right space (here `V`) we can use `CreateTimeRestrictedGF` which creates a new `GridFunction` and sets the values according to the space-time `GridFunction` restricted to a specific reference time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gfu = GridFunction(st_fes)\n",
    "u_last = CreateTimeRestrictedGF(gfu, 1)\n",
    "u, v = st_fes.TnT()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The exact solution is as stated above. With the help of a series of `GridFunction` we can visualize the function as a function in space for different values of `tref` that can be sampled by the `multidim`-component and animated:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u_exact = sin(pi * t) * sin(pi * x)**2 * sin(pi * y)**2\n",
    "from helper import ProjectOnMultiDimGF\n",
    "u_exact_to_show = ProjectOnMultiDimGF(sin(pi * tref) * sin(pi * x)**2 * sin(pi * y)**2,mesh,order=3,sampling=8)\n",
    "Draw(u_exact_to_show,mesh,\"u_exact\",autoscale=False,min=0,max=1, interpolate_multidim=True, animate=True, deformation=True)\n",
    "# The multidim-option now samples within the time interval [0,1]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The right-hand side $f$ is calculated accordingly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coeff_f = u_exact.Diff(t) - (u_exact.Diff(x).Diff(x) + u_exact.Diff(y).Diff(y))\n",
    "coeff_f = coeff_f.Compile()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variational formulation is derived from partial integration (in space) of the problem given above. We introduce the `DifferentialSymbol`s "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dxt = delta_t * dxtref(mesh, time_order=2)\n",
    "dxold = dmesh(mesh, tref=0)\n",
    "dxnew = dmesh(mesh, tref=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, `dxt` corresponds to $\\int_Q \\cdot d(x,t) = \\int_{t^{n-1}}^{t^n} \\int_{\\Omega} \\cdot d(x,t)$ whereas `dxold` and `dxnew` are spatial integrals on $\\Omega$ with however fixed reference time values 0 or 1 for the evaluation of space-time functions. With those, we can write the variational formulation as: Find $u \\in W_h = V_h \\times \\mathcal{P}^1([t^{n-1},t^n])$ so that\n",
    "$$\n",
    " \\int_{Q} (\\partial_t u v + \\nabla u \\cdot \\nabla v) ~~d(x,t) + \\int_{\\Omega} u_-(\\cdot,t^{n-1}) v_-(\\cdot,t^{n-1}) dx = \\int_Q f v ~ d(x,t) + \\int_{\\Omega} u_+(\\cdot,t^{n-1}) v_-(\\cdot,t^{n-1}) dx\\quad\\forall v \\in W_h\n",
    "$$\n",
    "where $u_+(\\cdot,t^{n-1})$ is the value from the last time step (initial value)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = BilinearForm(st_fes, symmetric=False)\n",
    "a += (dt(u) * v + grad(u) * grad(v)) * dxt\n",
    "a += u * v * dxold\n",
    "a.Assemble()\n",
    "ainv = a.mat.Inverse(st_fes.FreeDofs())\n",
    "\n",
    "f = LinearForm(st_fes)\n",
    "f += coeff_f * v * dxt\n",
    "f += u_last * v * dxold"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the l.h.s. does not depend on the time step. We hence assemble and factorize it only once.\n",
    "We use `u_last` as the `GridFunction` to hold the initial values. To initialize it, we take the space-time function `u_exact` and restrict it to `tref=0` using `fix_tref`. Not that `u_exact` depends on space and time whereas `u_last` is merely a spatial function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u_last.Set(fix_tref(u_exact, 0))\n",
    "told.Set(0.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For visualization purposes we store the initial value (and later add the solutions after each time step):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gfut = GridFunction(u_last.space, multidim=0)\n",
    "gfut.AddMultiDimComponent(u_last.vec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, the problem is solved in a time-stepping manner, which includes:\n",
    " * Assemble of r.h.s.\n",
    " * solve time step solution \n",
    " * extract solution on new time level (`RestrictGFInTime` extracts a spatial `GridFunction` from a space-time one for a fixed value `tref)\n",
    " * storing time step solution for later visualization (see last cell)\n",
    " * measuring the error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scene = Draw(u_last, mesh,\"u\", autoscale=False,min=0,max=1,deformation=True)\n",
    "timestep = 0\n",
    "while tend - told.Get() > delta_t / 2:\n",
    "    if timestep % 4 == 0:\n",
    "        gfut.AddMultiDimComponent(u_last.vec)    \n",
    "    timestep += 1\n",
    "    f.Assemble()\n",
    "    gfu.vec.data = ainv * f.vec\n",
    "    RestrictGFInTime(spacetime_gf=gfu, reference_time=1.0, space_gf=u_last)\n",
    "    l2error = sqrt(Integrate((u_exact - gfu)**2 * dxnew, mesh))\n",
    "    scene.Redraw()\n",
    "    told.Set(told.Get() + delta_t)\n",
    "    print(\"\\rt = {0:12.9f}, L2 error = {1:12.9e}\".format(told.Get(), l2error), end=\"\")\n",
    "gfut.AddMultiDimComponent(u_last.vec)    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we can visualize the time evoluation again. Note also the controls under `Multidim` (animate, speed, t)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw (gfut, mesh, interpolate_multidim=True, animate=True, deformation=True, min=0, max=1, autoscale=False)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.1-final"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}