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

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.dev0

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 0x7faa690a66b0>

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].

/usr/lib/python3/dist-packages/scipy/__init__.py:146: UserWarning: A NumPy version >=1.17.3 and <1.25.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"

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 += (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)
timestep = 0
while tend - told.Get() > delta_t / 2:
    if timestep % 4 == 0:
        gfut.AddMultiDimComponent(u_last.vec)
    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="")
gfut.AddMultiDimComponent(u_last.vec)

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