II. Maps, Tents & Nonlinear Equations
Jay Gopalakrishnan (gjay@pdx.edu)
July 2023 NGSolve User Meeting Lecture: Part II
A drawback of a spacetime tent region is that it is not a tensor product of a time interval with a spatial domain. This may seem to make timestepping algorithms and well-known spatial discretizations inapplicable on tents. But MTP (Mapped Tent Pitching) schemes overcame this drawback by mapping tents to tensor product domains. It allows the use of standard discontinuous Galerkin (DG) techniques after mapping tents. Once the mapped equation is understood,
you can quickly implement your favorite high-order spatial DG discretization in ngstents
while automatically leveraging the local time stepping capabilities of unstructured tent subdivisions.
Maps
The above-mentioned tent mapping is illustrated below. A tent \(\newcommand{\vtx}{{\mathtt{v}}} K_{\vtx}\) pitched atop a mesh vertex \(\vtx\) is shown as the image of a map \(\Phi\) from a tensor product spacetime cylinder \(\hat{K}_{\vtx}\). In the cylinder, instead of the physical time \(t\), we use a pseudotime \(\tau\) variable lying in the unit interval \([0, 1]\).
The map \(\Phi\) preserves space, but changes time per the next formula.
where \(\varphi\) is given as follows. Suppose the tent top and bottom are graphs of (continuous piecewise linear) functions \(\vphi_{\text{top}}(x)\) and \(\vphi_{\text{bot}}(x)\), respectively. Then
Let \(\delta(x) = \vphi_{\text{top}}(x) - \vphi_{\text{bot}}(x).\) Suppose we are given a hyperbolic system of \(L\) equations \(\partial_t U + \divx f(U ) = 0\) for a function \(U\) with \(L\) components, where the divergence is taken row wise on the \(L \times N\) matrix \(f(U)\) in \(N\)-spatial dimensions. The following result, proved in a [2017 paper], shapes the ensuing implementation details.
Transformation of the hyperbolic system:
A function \(U : K_{\vtx} \to \mathbb{R}^L\) satisfies
if and only if \(\; u = U \circ \Phi: \hat K_{\vtx} \to \mathbb{R}^L\) satisfies
Since the second equation is on the tensor product cylinder \(\hat{K}_{\vtx}\), we can use standard spatial discretizations there. The ngstents
package allows you to specify what DG method you want to apply, as we shall see next. Other details (such as mapping the solution back to the physical tent from the cylinder, tent dependencies and asynchronous parallel tent marching) are taken care of internally so you can focus on the discretization. In other words, ngstents
enables you to
experiment with very general hyperbolic systems and DG schemes by letting you provide your own \(f\) and your own DG schemes.
Inputs
The first required input is a function that returns the hyperbolic flux \(f(u) \in \mathbb{R}^{L \times N}\), given a state vector \(u \in \mathbb R^L\).
def HyperbolicFlux(u):
"""
Given L-vector u, compute and return L x N matrix f(u).
"""
Next, consider DG discretizations of the following form obtained by multiplying the mapped equation by a DG test function \(v\) and integrating by parts. Namely, a DG approximation \(u_h\) of \(u\) solves
where the “numerical flux” is that of your favorite DG method. That last term can be rewritten using facets. Each interior facet \(F \in \mathcal{F}_{int}\) is oriented by a unit normal \(n_F\). As one crosses \(F\) in the \(n_F\) direction, the value of a DG function \(u\) jumps from \(u^-\) to \(u^+\), i.e., \(u^\pm(x) = \lim_{s\to 0, s>0} u(x \pm s n_F)\) for \(x \in F\). On facets contained in the global boundary, \(n_F\) is always the outward unit normal \(n\). Using this notation,
where
and the volume term is \(T_{vol} = \sum_{K \subset \ov} \int_K \delta f(\uh_h) : \gradx v.\) Users can define their preferred DG method by providing definitions of \(\hat{f}(n_F, u^+, u^-)\) and \(\hat f_{bdr}( n, u^{bdr}, u)\) in python functions like those below.
def NumericalFlux(u_minus, u_plus):
"""
Given uminus and uplus on an interior facet F, which are
respectively traces from F's adjacent elements
from which and to which the facet's normal nF points,
compute and return fhat(n, u_minus, u_plus).
"""
def BoundaryFlux(u):
"""
Given u on a boundary facet, compute fhat(n, ubdr, u)
using known boundary data ubdr. Return it as a dictionary
with key-value pairs like
boundary_name: value of fhat(n, ubdr, u) at this named boundary
"""
Next, we come to an input function that is quite specific to the tent-based nature of the schemes. The user is required to provide the inverse of the map
The reason for it, without getting too lost in the weeds, is that it is needed to transform the equation, as derived previously,
into the more standard timesteppable form
def InverseMap(y):
"""
Given y, solve the equation y = u - f(u) grad(phi)
for u and return u.
"""
The quantity \(\gradx(\vphi)\) needed to solve the above equation is available as a data attribute of TentSlab
. The facet and boundary unit normal \(n\), usually needed in numerical flux calculations is provided by ngsolve
. These are usually enough to provide implementations of the above functions. These functions are then provided as input to create an MTP solver for linear hyperbolic systems. For nonlinear hyperbolic equations, a few more bells and whistles are needed for
stable shock capturing, as we will see later. But first, consider the following simple linear example.
Linear advection
Given a divergence-free vector field \(b\), let us solve the linear advection equation
using tents. We proceed to set up the above-mentioned four input functions for this example. They will be provided to an object of class ConservationLaw
, which we now import:
[1]:
from ngstents.conslaw import ConservationLaw
We repeat same workflow to pitch tents as in the earlier tutorial:
[2]:
from ngstents import TentSlab
import ngsolve as ng
from ngsolve import CF, L2, sqrt, exp, x, y, IfPos, InnerProduct
from netgen.geom2d import unit_square
from ngsolve.webgui import Draw
[3]:
mesh = ng.Mesh(unit_square.GenerateMesh(maxh=0.05))
dt = 0.02
ts = TentSlab(mesh)
ts.SetMaxWavespeed(1)
ts.PitchTents(dt=dt, local_ct=True, global_ct=0.5);
Next, we set the advection vector field \(b\) and the initial function.
[4]:
b = CF((1, 0.2), dims=(1, 2))
u0 = CF(exp(-50 * ((x - 0.3)**2 + (y - 0.3)**2)))
Draw(u0, mesh);
Now we are ready to define the four required input functions for solving this problem using tents.
[5]:
n = ng.specialcf.normal(mesh.dim)
def HyperbolicFlux_Advection(u): # f(u) = b * u
return CF(b * u, dims=(1, 2))
def NumericalFlux_Advection(um, up): # set upwind flux
bn = b * n
return IfPos(bn, bn * um, bn * up)
def BoundaryFlux_Advection(u):
return mesh.BoundaryCF({
"left": 0, # inflow
"bottom": 0, # inflow
"right": b * n * u, # outflow
"top": b * n * u # outflow
})
def InverseMap_Advection(y):
# solve for u from the linear eq y = u - u * b * grad(phi)
return y / (1 - InnerProduct(b, ts.gradphi))
Providing these inputs to a conservation law object, and using its propagation method, we obtain the solution to the advection equation.
[6]:
V = L2(mesh, order=4)
u = ng.GridFunction(V)
cl = ConservationLaw(u, ts,
flux=HyperbolicFlux_Advection,
numflux=NumericalFlux_Advection,
inversemap=InverseMap_Advection)
cl.SetInitial(u0)
cl.SetBoundaryCF(BoundaryFlux_Advection(cl.u_minus))
cl.SetTentSolver(substeps=10)
t = 0; tend = 0.6
scene = Draw(u)
with ng.TaskManager():
while t < tend - dt / 2:
cl.Propagate()
t += dt
scene.Redraw()
set up 2-stage (second order) SARK timestepping with 10 substeps/tent
Two-dimensional Burger’s equation
Let us now turn to the more complex case of a nonlinear hyperbolic equation. We illustrate the usage of MTP techniques considering the example of a two-dimensional Burger’s equation,
In addition to the above-mentioned four required input functions, we may now specify four more input functions that allow for stable simulations even in the presence of shocks. The literature contains many shock capturing techniques (like flux and slope limiters) for nonlinear conservation laws. We implemented one called the entropy viscosity regularization of Guermond, Pasquetti and Popov’s [2011 paper]. They suggest modifying numerical schemes by adding small amounts of artificial viscosity in selected locations to avoid spurious oscillations near shocks.
How are these locations selected? Recall that a scalar function \(\Ec(u)\) is called an entropy if there exists a vector field called entropy flux \(\Fc(u)\) such that any smooth solution \(u\) of the hyperbolic system satisfies the entropy conservation law
For a nonsmooth solution \(u\), only the entropy admissibility inequality
can be expected to hold rather than the equality. By computing an entropy residual \(r_h \approx \partial_t \Ec(u)_h + \divx \Fc(u_h)\) from a numerical solution \(u_h\), regions of entropy production are identified as locations for viscosity addition. The following four functions are required from the user to turn on entropy viscosity regularization.
def Entropy(u):
"""
Compute and return a known entropy function E(u).
"""
def EntropyFlux(u):
"""
Compute the entropy flux F(u) for the above entropy E(u).
"""
def NumericalEntropyFlux(u_minus, u_plus):
"""
Return a DG numerical flux of a discretization of the entropy
conservation law to be used for entropy residual calculation.
Used only on interior facets.
"""
def ViscosityCoefficient(u, res):
"""
Given an approximate solution "u" and its entropy residual "res",
return the artificial viscosity coefficient to be added.
"""
For the Burger’s equation, a well known entropy pair is
This together with upwind DG fluxes are set in the following functions.
[7]:
# PITCH TENTS
mesh = ng.Mesh(unit_square.GenerateMesh(maxh=0.05))
dt = 0.01
ts = TentSlab(mesh) # , method="edge")
ts.SetMaxWavespeed(2)
ts.PitchTents(dt=dt, local_ct=True, global_ct=0.5);
n = ng.specialcf.normal(mesh.dim)
h = ng.specialcf.mesh_size
[8]:
# SET 8 INPUT FUNCTIONS FOR 2D BURGER'S EQUATION
def HyperbolicFlux_Burgers(u):
return CF((u**2 / 2, u**2 / 2), dims=(1, 2))
def NumericalFlux_Burgers(u_minus, u_plus):
uhat = 0.5 * (u_minus + u_plus)
dfn = uhat * (n[0] + n[1]) # dfn = f'(u) n
# Set upwind flux determined by sgn(f'(u) n)
return IfPos(dfn,
HyperbolicFlux_Burgers(u_minus) * n,
HyperbolicFlux_Burgers(u_plus) * n)
def BoundaryFlux_Burgers(u):
bdrycf = {bdr: NumericalFlux_Burgers(u, 0) for bdr in mesh.GetBoundaries()}
return mesh.BoundaryCF(bdrycf)
def InverseMap_Burgers(y):
# Solve the quadratic equation
# y = u - (1/2)[u^2, u^2] grad(phi), or equivalently,
# y = u - u^2 * d, where d is as follows:
d = ts.gradphi[0] + ts.gradphi[1]
return 2 * y / (1 + sqrt(1 - 2 * d * y))
def Entropy_Burgers(u):
return u**2 / 2
def EntropyFlux_Burgers(u):
return CF((u**3 / 3, u**3 / 3), dims=(1, 2))
def NumericalEntropyFlux_Burgers(u_minus, u_plus):
uhat = 0.5 * (u_minus + u_plus)
dFn = uhat**2 * (n[0] + n[1]) # F'(u) * n
# Set upwind flux determined by sgn(F'(u) n)
return IfPos(dFn,
EntropyFlux_Burgers(u_minus) * n,
EntropyFlux_Burgers(u_plus) * n)
def ViscosityCoefficient_Burgers(u, res):
R = IfPos(res, res, -res)
Eu = Entropy_Burgers(u)
E = IfPos(Eu - 1e-13, Eu, 1e-13)
entropy_viscosity = (h / p)**2 * IfPos(R - 1e-13, R/E, 0)
entropy_viscosmax = (h / p) * IfPos(u, u, -u)
return IfPos(entropy_viscosmax - entropy_viscosity, entropy_viscosity,
entropy_viscosmax)
[9]:
# INITIALIZE CONSERVATION LAW OBJECT
p = 4
V = ng.L2(mesh, order=p)
u = ng.GridFunction(V)
u0 = CF(exp(-50 * ((x - 0.3)**2 + (y - 0.3)**2)))
burg = ConservationLaw(u, ts,
flux=HyperbolicFlux_Burgers,
numflux=NumericalFlux_Burgers,
inversemap=InverseMap_Burgers,
entropy=Entropy_Burgers,
entropyflux=EntropyFlux_Burgers,
numentropyflux=NumericalEntropyFlux_Burgers,
visccoeff=ViscosityCoefficient_Burgers)
burg.SetInitial(u0)
burg.SetBoundaryCF(BoundaryFlux_Burgers(burg.u_minus))
burg.SetTentSolver(substeps=10)
set up 2-stage (second order) SARK timestepping with 10 substeps/tent
[10]:
ut = ng.GridFunction(V, multidim=0)
ut.AddMultiDimComponent(burg.sol.vec)
entropyresidual = ng.GridFunction(V, multidim=0)
entropyresidual.AddMultiDimComponent(burg.res.vec)
At this point we have set up all the ingredients for solving the nonlinear Burger’s equation using tents. Both the solution and the entropy residual are stored in a multidimensional grid function for visualization after the computations.
[11]:
t = 0
tend = 0.6
scene = Draw(burg.sol)
with ng.TaskManager():
while t < tend - dt / 2:
burg.Propagate()
t += dt
scene.Redraw()
ut.AddMultiDimComponent(burg.sol.vec)
entropyresidual.AddMultiDimComponent(burg.res.vec)
[12]:
Draw(ut, mesh, autoscale=False, min=0, max=0.7, interpolate_multidim=True, animate=True);
Here is an animation of the computed entropy residual. It is notably large in regions of sharp transitions.
[13]:
Draw(entropyresidual, mesh, interpolate_multidim=True, animate=True,
autoscale=False, min=-1e-6, max=1e-6);
Conclusion
Here is a summary of what we have seen above:
Mapping spacetime tents to tensor product domains
Using DG spatial discretizations on the tensor product domain
4 python input functions for spatial discretization of linear problems
8 python input functions for solving nonlinear hyperbolic equations
Solving Burger’s equation with entropy viscosity regularization