II. Maps, Tents & Nonlinear Equations

Jay Gopalakrishnan (gjay@pdx.edu)

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.

\[\begin{split}\begin{align*} \newcommand{\d}{\partial} \newcommand{\divxh}{\hat{\mathrm{div}}_x} \newcommand{\divhx}{\text{div}_{x}} \newcommand{\divx}{\text{div}_x} \newcommand{\gradx}{\text{grad}_x} \newcommand{\gradhx}{\mathop{\text{grad}_{x}}} \newcommand{\vphi}{\varphi} \Phi \begin{pmatrix} {{x}} \\ \,{\tau} \end{pmatrix} & = \begin{pmatrix} {{x}} \\ \varphi( {{x}}, {\tau}) \end{pmatrix}, \end{align*}\end{split}\]

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

\[ \vphi(x, \tau) = (1-{\tau}) \vphi_{\text{bot}} + {\tau} \vphi_{\text{top}}. \]

This change of variable will change the equation we want to solve. Suppose we are initially given a hyperbolic system of \(L\) equations \( \newcommand{\vtx}{\text{v}} \newcommand{\d}{\partial} \newcommand{\divxh}{\hat{\mathrm{div}}_x} \newcommand{\divhx}{\text{div}_{x}} \newcommand{\divx}{\text{div}_x} \newcommand{\gradx}{\text{grad}_x} \newcommand{\vphi}{\varphi} \newcommand{\gradhx}{\mathop{\text{grad}_{x}}} \)

\[ \partial_t U + \divx f(U ) = 0 \]

for some function \(U\) with \(L\) components on the tent \(K_{\vtx}\). Here the divergence is taken row wise on the \(L \times N\) matrix \(f(U)\) in \(N\)-spatial dimensions. The transformed variable is

\[ u = U \circ \Phi: \hat K_{\vtx} \to \mathbb{R}^L. \]

The following result, proved in a [2017 paper], shows how the equation transforms under the variable change, and shapes the ensuing implementation details.

Transformation of the hyperbolic system: The equation

\[ \def\d{\partial} \newcommand{\vtx}{\text{v}} \newcommand{\d}{\partial} \newcommand{\divxh}{\hat{\mathrm{div}}_x} \newcommand{\divhx}{\text{div}_{x}} \newcommand{\divx}{\text{div}_x} \newcommand{\gradx}{\text{grad}_x} \newcommand{\vphi}{\varphi} \newcommand{\gradhx}{\mathop{\text{grad}_{x}}} \frac{\d U}{\d t} + \divx f( {U}) =0,\qquad \text { on } K_{\vtx}, \]

holds if and only if \(\; u = U \circ \Phi\) on \(\hat K_{\vtx}\) satisfies

\[ \frac{\d }{\d {\tau}} \big[ u - f(u)\, \gradx \vphi \big] + \divhx {\delta} f(u ) = 0, \qquad \text { on } \hat K_{\vtx}. \]

where \(\delta(x) = \vphi_{\text{top}}(x) - \vphi_{\text{bot}}(x).\)

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

\[\begin{split}\begin{align*} \newcommand{\uh}{{{u}}} \newcommand{\ov}{{\Omega_{\vtx}}} \newcommand{\oh}{{\varOmega_h}} % \int_{\Omega_{\vtx}} & \frac{\partial }{\partial {\tau}} \big[ \uh_h - f(\uh_h) \gradx \varphi\big] \cdot v = \sum_{K \subset \ov} \bigg[ \int_K \delta f(\uh_h) : \gradx v - \int_{\partial K} \left(\begin{smallmatrix} \text{numerical}\\ \text{flux} \end{smallmatrix}\right) \cdot \,v \,\delta\,\bigg], \end{align*}\end{split}\]

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,

\[\newcommand{\uh}{{{u}}} \newcommand{\ov}{{\Omega_{\vtx}}} \newcommand{\oh}{{\varOmega_h}} \frac{\partial }{\partial {\tau}} \int_{\Omega_{\vtx}} \big[ \uh_h - f(\uh_h) \gradx \varphi\big] \cdot v =T_{vol} + T_{int} + T_{bdr} \]

where

\[ T_{int} = -\sum_{F\in \mathcal{F}_{int}} \int_F \delta \;\hat{f}(n_F, u_h^+, u_h^-) \;(v^+ - v^-), \qquad T_{bdr} = \sum_{F \in \mathcal{F}_{bdr}} \int_F \delta \; \hat f_{bdr}( n, u^{bdr}, u_h)\; v \]

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

\[ u_h \longmapsto y_h = u_h - f(u_h) \gradx \vphi. \]

The reason for it, without getting too lost in the weeds, is that it is needed to transform the equation, as derived previously,

\[ \frac{d }{d\tau} \int_{\Omega_{\vtx}} \big[ \uh_h - f(\uh_h) \gradx \varphi\big] \cdot v = \ldots {\text{some function of $u_h$}\ldots } \]

into the more standard timesteppable form

\[ \frac{d }{d {\tau}} \int_{\Omega_{\vtx}} y_h \cdot v = \ldots \text{some function of } y_h \ldots. \]

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

\[ \partial_t u - \divx( b u ) = 0 \]

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:

from ngstents.conslaw import ConservationLaw

We repeat previous workflow to pitch tents as in the earlier tutorial:

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

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.

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.

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.

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,

\[\newcommand{\Ec}{\mathcal{E}} \newcommand{\Fc}{\mathcal{F}} \partial_t U + \divx \left(\frac 1 2 [U^2, U^2]\right) = 0. \]

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 \(\newcommand{\Ec}{\mathcal{E}} \newcommand{\Fc}{\mathcal{F}}\) \(\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

\[ \newcommand{\Ec}{\mathcal{E}} \newcommand{\Fc}{\mathcal{F}} \partial_t \Ec(u) + \divx \Fc(u) = 0. \]

For a nonsmooth solution \(u\), only the entropy admissibility inequality

\[ \partial_t \Ec(u) + \divx \Fc(u) \le 0 \]

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 2D Burger’s equation, a well known entropy pair is

\[ \Ec(u) = \frac 1 2 u^2, \qquad \Fc(u) = \frac 1 3 [u^3, u^3]. \]

This together with upwind DG fluxes are set in the following functions.

# 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

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

# 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

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.

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)

Draw(ut, mesh, autoscale=False, min=0, max=0.7, interpolate_multidim=True, animate=True, scale=0.2, deformation=True, euler_angles=[-35,-1, -93]);

Here is an animation of the computed entropy residual. It is notably large in regions of sharp transitions.

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 2D Burger’s equation with entropy viscosity regularization