This page was generated from unit-8.4-spacetime_fitted/spacetime_fitted.ipynb.

In this tutorial we make use of ipywidgets. You can (in the ideal case) simply install them with executing the next cell. You will however need to restart jupyter afterwards (unless you have it installed already).

!pip3 install --user ipywidgets
!jupyter nbextension enable --py widgetsnbextension

Requirement already satisfied: ipywidgets in /usr/local/lib/python3.10/dist-packages (8.0.6)
Requirement already satisfied: traitlets>=4.3.1 in /usr/local/lib/python3.10/dist-packages (from ipywidgets) (5.9.0)
Requirement already satisfied: widgetsnbextension~=4.0.7 in /usr/local/lib/python3.10/dist-packages (from ipywidgets) (4.0.7)
Requirement already satisfied: ipython>=6.1.0 in /usr/local/lib/python3.10/dist-packages (from ipywidgets) (8.14.0)
Requirement already satisfied: jupyterlab-widgets~=3.0.7 in /usr/local/lib/python3.10/dist-packages (from ipywidgets) (3.0.7)
Requirement already satisfied: ipykernel>=4.5.1 in /usr/local/lib/python3.10/dist-packages (from ipywidgets) (6.23.3)
Requirement already satisfied: comm>=0.1.1 in /usr/local/lib/python3.10/dist-packages (from ipykernel>=4.5.1->ipywidgets) (0.1.3)
Requirement already satisfied: pyzmq>=20 in /usr/local/lib/python3.10/dist-packages (from ipykernel>=4.5.1->ipywidgets) (25.1.0)
Requirement already satisfied: tornado>=6.1 in /usr/local/lib/python3.10/dist-packages (from ipykernel>=4.5.1->ipywidgets) (6.3.2)
Requirement already satisfied: matplotlib-inline>=0.1 in /usr/local/lib/python3.10/dist-packages (from ipykernel>=4.5.1->ipywidgets) (0.1.6)
Requirement already satisfied: jupyter-core!=5.0.*,>=4.12 in /usr/local/lib/python3.10/dist-packages (from ipykernel>=4.5.1->ipywidgets) (5.3.1)
Requirement already satisfied: debugpy>=1.6.5 in /usr/local/lib/python3.10/dist-packages (from ipykernel>=4.5.1->ipywidgets) (1.6.7)
Requirement already satisfied: packaging in /usr/lib/python3/dist-packages (from ipykernel>=4.5.1->ipywidgets) (21.3)
Requirement already satisfied: jupyter-client>=6.1.12 in /usr/local/lib/python3.10/dist-packages (from ipykernel>=4.5.1->ipywidgets) (8.3.0)
Requirement already satisfied: psutil in /usr/local/lib/python3.10/dist-packages (from ipykernel>=4.5.1->ipywidgets) (5.9.5)
Requirement already satisfied: nest-asyncio in /usr/local/lib/python3.10/dist-packages (from ipykernel>=4.5.1->ipywidgets) (1.5.6)
Requirement already satisfied: pygments>=2.4.0 in /usr/lib/python3/dist-packages (from ipython>=6.1.0->ipywidgets) (2.11.2)
Requirement already satisfied: jedi>=0.16 in /usr/local/lib/python3.10/dist-packages (from ipython>=6.1.0->ipywidgets) (0.18.2)
Requirement already satisfied: pexpect>4.3 in /usr/local/lib/python3.10/dist-packages (from ipython>=6.1.0->ipywidgets) (4.8.0)
Requirement already satisfied: decorator in /usr/lib/python3/dist-packages (from ipython>=6.1.0->ipywidgets) (4.4.2)
Requirement already satisfied: stack-data in /usr/local/lib/python3.10/dist-packages (from ipython>=6.1.0->ipywidgets) (0.6.2)
Requirement already satisfied: pickleshare in /usr/local/lib/python3.10/dist-packages (from ipython>=6.1.0->ipywidgets) (0.7.5)
Requirement already satisfied: backcall in /usr/local/lib/python3.10/dist-packages (from ipython>=6.1.0->ipywidgets) (0.2.0)
Requirement already satisfied: prompt-toolkit!=3.0.37,<3.1.0,>=3.0.30 in /usr/local/lib/python3.10/dist-packages (from ipython>=6.1.0->ipywidgets) (3.0.38)
Requirement already satisfied: parso<0.9.0,>=0.8.0 in /usr/local/lib/python3.10/dist-packages (from jedi>=0.16->ipython>=6.1.0->ipywidgets) (0.8.3)
Requirement already satisfied: python-dateutil>=2.8.2 in /usr/local/lib/python3.10/dist-packages (from jupyter-client>=6.1.12->ipykernel>=4.5.1->ipywidgets) (2.8.2)
Requirement already satisfied: platformdirs>=2.5 in /usr/local/lib/python3.10/dist-packages (from jupyter-core!=5.0.*,>=4.12->ipykernel>=4.5.1->ipywidgets) (3.8.0)
Requirement already satisfied: ptyprocess>=0.5 in /usr/local/lib/python3.10/dist-packages (from pexpect>4.3->ipython>=6.1.0->ipywidgets) (0.7.0)
Requirement already satisfied: wcwidth in /usr/local/lib/python3.10/dist-packages (from prompt-toolkit!=3.0.37,<3.1.0,>=3.0.30->ipython>=6.1.0->ipywidgets) (0.2.6)
Requirement already satisfied: asttokens>=2.1.0 in /usr/local/lib/python3.10/dist-packages (from stack-data->ipython>=6.1.0->ipywidgets) (2.2.1)
Requirement already satisfied: pure-eval in /usr/local/lib/python3.10/dist-packages (from stack-data->ipython>=6.1.0->ipywidgets) (0.2.2)
Requirement already satisfied: executing>=1.2.0 in /usr/local/lib/python3.10/dist-packages (from stack-data->ipython>=6.1.0->ipywidgets) (1.2.0)
Requirement already satisfied: six in /usr/lib/python3/dist-packages (from asttokens>=2.1.0->stack-data->ipython>=6.1.0->ipywidgets) (1.16.0)
WARNING: Running pip as the 'root' user can result in broken permissions and conflicting behaviour with the system package manager. It is recommended to use a virtual environment instead: https://pip.pypa.io/warnings/venv
Enabling notebook extension jupyter-js-widgets/extension...
      - Validating: OK

8.4 Space-time discretizations on fitted geometry

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:

\[\begin{split}\left\{ \begin{aligned} \partial_t u - \Delta u &= f \quad \text{ in } \Omega \text{ for all } t \in [0,1], & \\ ~ \partial_{\mathbf{n}} u &= 0 \quad \text{ on } \partial \Omega, & \\ u &= u_0 \quad \text{at } t=0. & \\ \end{aligned}\right.\end{split}\]

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\)

from netgen.geom2d import unit_square
from ngsolve import *
from xfem import *
import time
from math import pi
ngsglobals.msg_level = 1

importing ngsxfem-2.1.2301

We opt for a first-order method in space and time and choose an appropriate timestep.

# Space finite element order
order = 1
# Time finite element order
k_t = 1
# Final simulation time
tend = 1.0
# Time step
delta_t = 1 / 32

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:

def dt(u):
    return 1.0 / delta_t * dtref(u)

tnew = 0
told = Parameter(0)
t = told + delta_t * tref

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.

mesh = Mesh(unit_square.GenerateMesh(maxh=0.05, quad_dominated=False))

V = H1(mesh, order=order, dirichlet=".*")
tfe = ScalarTimeFE(k_t)
st_fes = tfe * V # tensor product of a space finite element space and a 1D finite element

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.

gfu = GridFunction(st_fes)
u_last = CreateTimeRestrictedGF(gfu, 1)
u, v = st_fes.TnT()

The exact solution is as stated above. With the TimeSlider_Draw we can visualize the function as a function in space for different values of tref

u_exact = sin(pi * t) * sin(pi * x)**2 * sin(pi * y)**2
TimeSlider_Draw(sin(pi * tref) * sin(pi * x)**2 * sin(pi * y)**2,mesh,"u_exact", autoscale=False,min=0,max=1)
#In order to draw the function u_exact for t in [0,1] instead of t in [0, delta t], we rewrite the function

<function xfem.TimeSlider_Draw.<locals>.UpdateTime(time)>

The right-hand side \(f\) is calculated accordingly and can also be visualized similarly:

coeff_f = u_exact.Diff(t) - (u_exact.Diff(x).Diff(x) + u_exact.Diff(y).Diff(y))
coeff_f = coeff_f.Compile()
TimeSlider_Draw(coeff_f,mesh,"coeff_f", autoscale=False,min=-2,max=8)

Warning: differentiationg by a variable not marked as Variable,
might be optimized out. Call MakeVariable for differentiation CF

<function xfem.TimeSlider_Draw.<locals>.UpdateTime(time)>

The variational formulation is derived from partial integration (in space) of the problem given above. We introduce the DifferentialSymbols

dxt = delta_t * dxtref(mesh, time_order=2)
dxold = dmesh(mesh, tref=0)
dxnew = dmesh(mesh, tref=1)

Here, dxt corresponds to $ \intQ :nbsphinx-math:`cdot `d(x,t) = :nbsphinx-math:`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

\[\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\]

where \(u_+(\cdot,t^{n-1})\) is the value from the last time step (initial value).

a = BilinearForm(st_fes, symmetric=False)
a += (dt(u) * v + grad(u) * grad(v)) * dxt
a += u * v * dxold
a.Assemble()
ainv = a.mat.Inverse(st_fes.FreeDofs())

f = LinearForm(st_fes)
f += coeff_f * v * dxt
f += u_last * v * dxold

Note that the l.h.s. does not depend on the time step. We hence assemble and factorize it only once. 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.

u_last.Set(fix_tref(u_exact, 0))
told.Set(0.)

For visualization purposes we store the initial value (and later add the solutions after each time step):

gfut = GridFunction(u_last.space, multidim=0)
gfut.AddMultiDimComponent(u_last.vec)

Finally, the problem is solved in a time-stepping manner, which includes: * Assemble of r.h.s. * solve time step solution * extract solution on new time level (RestrictGFInTime extracts a spatial GridFunction from a space-time one for a fixed value `tref) * storing time step solution for later visualization (see last cell) * measuring the error

scene = Draw(u_last, mesh,"u", autoscale=False,min=0,max=1,deformation=True)
while tend - told.Get() > delta_t / 2:
    f.Assemble()
    gfu.vec.data = ainv * f.vec
    RestrictGFInTime(spacetime_gf=gfu, reference_time=1.0, space_gf=u_last)
    gfut.AddMultiDimComponent(u_last.vec)
    l2error = sqrt(Integrate((u_exact - gfu)**2 * dxnew, mesh))
    scene.Redraw()
    told.Set(told.Get() + delta_t)
    print("\rt = {0:12.9f}, L2 error = {1:12.9e}".format(told.Get(), l2error), end="")

check type <class 'ngsolve.comp.GridFunction'>
check type <class 'ngsolve.comp.GridFunction'>

t =  1.000000000, L2 error = 1.840270280e-04

Next, we can visualize the time evoluation again. Note also the controls under Multidim (animate, speed, t).

Draw (gfut, mesh, interpolate_multidim=True, animate=True, deformation=True)

BaseWebGuiScene