# 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.2303


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
t.MakeVariable()

:

<ngsolve.fem.CoefficientFunction at 0x7f48d6701da0>


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 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:

:

u_exact = sin(pi * t) * sin(pi * x)**2 * sin(pi * y)**2
from helper import ProjectOnMultiDimGF
u_exact_to_show = ProjectOnMultiDimGF(sin(pi * tref) * sin(pi * x)**2 * sin(pi * y)**2,mesh,order=3,sampling=8)
Draw(u_exact_to_show,mesh,"u_exact",autoscale=False,min=0,max=1, interpolate_multidim=True, animate=True, deformation=True)
# The multidim-option now samples within the time interval [0,1].

:

BaseWebGuiScene


The right-hand side $$f$$ is calculated accordingly:

:

coeff_f = u_exact.Diff(t) - (u_exact.Diff(x).Diff(x) + u_exact.Diff(y).Diff(y))
coeff_f = coeff_f.Compile()


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

$\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 += 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)


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)
timestep = 0
while tend - told.Get() > delta_t / 2:
if timestep % 4 == 0:
timestep += 1
f.Assemble()
gfu.vec.data = ainv * f.vec
RestrictGFInTime(spacetime_gf=gfu, reference_time=1.0, space_gf=u_last)
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="")

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, min=0, max=1, autoscale=False)

:

BaseWebGuiScene
`