This page was generated from jupyter-files/unit-3.4-simplehyp/shallow2D.ipynb.

3.4 A Nonlinear conservation law: the shallow water equation in 2D

We consider the shallow water equations as an example of a nonlinear conservation law, i.e. we consider

\[\partial_t \mathbf{U} + \operatorname{div}(\mathbf{F} (\mathbf{U} )) = 0 \qquad in \qquad \Omega \times[0,T],\]

with

\[\mathbf{U} = (h, \mathbf{w}) = (h, hu, hv) = (\mathbf{u}_1, \mathbf{u}_2, \mathbf{u}_3)\]

and

\[\begin{split}\mathbf{F}(\mathbf{U}) = \left( \begin{array}{c@{}c@{\qquad}c@{}c} h u & & h v & \\ h u^2 &+ \frac12 g h^2 & h u v & \\ h u v & & h v^2 & + \frac12 g h^2 \end{array} \right) = \left( \begin{array}{cc} \mathbf{u}_2 & \mathbf{u}_3 \\ \frac{\mathbf{u}_2^2}{\mathbf{u}_1} + \frac12 g \mathbf{u}_1^2 & \frac{\mathbf{u}_2 \mathbf{u}_3}{\mathbf{u}_1} \\ \frac{\mathbf{u}_2 \mathbf{u}_3}{\mathbf{u}_1} & \frac{\mathbf{u}_3^2}{\mathbf{u}_1} + \frac12 g \mathbf{u}_1^2 \end{array} \right) = \left( \begin{array}{cc} h \mathbf{w}^T \\ h \mathbf{w} \otimes \mathbf{w} + \frac12 g h^2 \mathbf{I} \end{array} \right)\end{split}\]

Jacobian of the flux for shallow water:

\begin{align*} \mathbf{A}_1 & = \left( \begin{array}{ccc} 0 & 1 & 0 \\ - \frac{\mathbf{u}_2^2}{\mathbf{u}_1^2} + g \mathbf{u}_1 & 2 \frac{\mathbf{u}_2}{\mathbf{u}_1} & 0 \\ - \frac{\mathbf{u}_2 \mathbf{u}_3}{\mathbf{u}_1^2} & \frac{\mathbf{u}_3}{\mathbf{u}_1} & \frac{\mathbf{u}_2}{\mathbf{u}_1} \end{array} \right) = \left( \begin{array}{ccc} 0 & 1 & 0 \\ - u^2 + g h & 2 u & 0 \\ - u v & v & u \end{array} \right) \\ \mathbf{A}_2 & = \left( \begin{array}{ccc} 0 & 0 & 1 \\ - \frac{\mathbf{u}_2\mathbf{u}_3}{\mathbf{u}_1^2} & \frac{\mathbf{u}_3}{\mathbf{u}_1} & \frac{\mathbf{u}_2}{\mathbf{u}_1} \\ - \frac{\mathbf{u}_3^2}{\mathbf{u}_1^2} + g \mathbf{u}_1 & 0 & 2\frac{\mathbf{u}_3}{\mathbf{u}_1} \end{array} \right) = \left( \begin{array}{ccc} 0 & 0 & 1 \\ - uv & v & u \\ - v^2 + gh & 0 & 2 v \end{array} \right) \\ \mathbf{A}(\mathbf{u}, \mathbf{n}) & = n_1 \mathbf{A}_1 + n_2 \mathbf{A}_2 = \left( \begin{array}{ccc} 0 & n_1 & n_2 \\ - u \alpha - g h n_1 & \alpha + u n_1 & u n_2 \\ - v \alpha - g h n_2 & v n_1 & \alpha + v n_2 \end{array} \right), \quad \text{ with } \alpha = (\mathbf{u},\mathbf{n}), \end{align*}

spectrum:

\[\rho( \mathbf{A}(\mathbf{u}, \mathbf{n}) ) = \{ \alpha + \sqrt{g h}, \alpha, \alpha - \sqrt{g h} \}\]

from ngsolve import * # sloppy
import netgen.gui
%gui tk

The dam break problem (geometry)

from netgen.geom2d import SplineGeometry
geo = SplineGeometry()

pnts =[ (-12,-5), (-7,-5), (-5,-5), (-3,-5), (-3,-3),
        (-3,-1), (-1,-1), ( 0,-1), ( 1,-1), ( 3,-1),
        ( 3,-3), ( 3,-5), ( 5,-5), ( 7,-5), ( 12,-5),
        ( 12, 5), ( 7, 5), ( 5, 5), ( 3, 5), ( 3, 3),
        ( 3, 1), ( 1, 1), ( 0, 1), (-1, 1), (-3, 1),
        (-3, 3), (-3, 5), (-5, 5), (-7, 5), (-12, 5)]
pnts = [geo.AppendPoint(*pnt) for pnt in pnts]

curves = [[["line",0,1],"wall",1, 0], [["line",1,2],"wall",1, 0],
          [["spline3",2,3,4],"wall",1, 0],[["spline3",4,5,6],"wall",1, 0],[["line",6,7],"wall",1, 0],
          [["line",7,22],"dam",1, 2],     # <--- dam interface
          [["line",7,8],"wall",2, 0],[["spline3",8,9,10],"wall",2, 0],
          [["spline3",10,11,12],"wall",2, 0],[["line",12,13],"wall",2, 0],[["line",13,14],"wall",2, 0],
          [["line",14,15],"wall",2, 0], # <--- right boundary
          [["line",15,16],"wall",2, 0],[["line",16,17],"wall",2, 0],
          [["spline3",17,18,19],"wall",2, 0],[["spline3",19,20,21],"wall",2, 0],
          [["line",21,22],"wall",2, 0],[["line",22,23],"wall",1, 0],
          [["spline3",23,24,25],"wall",1, 0],[["spline3",25,26,27],"wall",1, 0],
          [["line",27,28],"wall",1, 0],[["line",28,29],"wall",1, 0],
          [["line",29,0],"wall",1, 0]]   # <--- left boundary

for c,bc,l,r in curves:
    geo.Append(c,bc=bc,leftdomain=l, rightdomain=r)
geo.SetMaterial(1,"upperlevel")
geo.SetMaterial(2,"lowerlevel")

mesh = Mesh(geo.GenerateMesh(maxh=2))
mesh.Curve(3)
Draw(mesh)

Approximation space

order = 2
fes = L2(mesh,order=order,dim=3)

initial and boundary conditions

U,V = fes.TnT() # "Trial" and "Test" function
h, hu, hv = U

# initial conditions
h0mat = {"upperlevel" : 10, "lowerlevel" : 2}
U0 = CoefficientFunction((CoefficientFunction([h0mat[mat] for mat in mesh.GetMaterials()]),0,0))

# boundary conditions
hbndreg = CoefficientFunction([{"wall" : h, "dam" : 0}[rg] for rg in mesh.GetBoundaries()])
hubndreg = CoefficientFunction([{"wall" : -hu, "dam" : 0}[rg] for rg in mesh.GetBoundaries()])
hvbndreg = CoefficientFunction([{"wall" : -hv, "dam" : 0}[rg] for rg in mesh.GetBoundaries()])

Ubnd = CoefficientFunction((hbndreg,hubndreg,hvbndreg))

# constant for gravitational force
g=1

Flux definition and numerical flux choice (Lax-Friedrich)

\[\begin{split}\mathbf{F}(\mathbf{U}) = \left( \begin{array}{c@{}c@{\qquad}c@{}c} h u & & h v & \\ h u^2 &+ \frac12 g h^2 & h u v & \\ h u v & & h v^2 & + \frac12 g h^2 \end{array} \right) = \left( \begin{array}{cc} \mathbf{u}_2 & \mathbf{u}_3 \\ \frac{\mathbf{u}_2^2}{\mathbf{u}_1} + \frac12 g \mathbf{u}_1^2 & \frac{\mathbf{u}_2 \mathbf{u}_3}{\mathbf{u}_1} \\ \frac{\mathbf{u}_2 \mathbf{u}_3}{\mathbf{u}_1} & \frac{\mathbf{u}_3^2}{\mathbf{u}_1} + \frac12 g \mathbf{u}_1^2 \end{array} \right) = \left( \begin{array}{cc} h \mathbf{w}^T \\ h \mathbf{w} \otimes \mathbf{w} + \frac12 g h^2 \mathbf{I} \end{array} \right)\end{split}\]

def F(U):
    h, hvx, hvy = U
    vx = hvx/h
    vy = hvy/h
    return CoefficientFunction(((hvx,hvy),
                                (hvx*vx + 0.5*g*h**2, hvx*vy),
                                (hvx*vy, hvy*vy + 0.5*g*h**2)),dims=(3,2))

\[\hat{\mathbf{F}}(\mathbf{U}) \cdot \mathbf{n} = \frac{(\mathbf{F}(\mathbf{U}^l)+\mathbf{F}(\mathbf{U}^r)}{2} \cdot \mathbf{n} + \frac{|\lambda|}{2} (\mathbf{U}^l - \mathbf{U}^r), \quad (\cdot^l: \text{current element}, \quad \cdot^r: \text{neighbor element})\]

n = specialcf.normal(mesh.dim)

def Max(u,v):
    return IfPos(u-v,u,v)
def Fmax(A,B):
    ha, hua, hva = A
    hb, hub, hvb = B
    vnorma = sqrt(hua**2+hva**2)/ha
    vnormb = sqrt(hub**2+hvb**2)/hb
    return Max(vnorma+sqrt(g*A[0]),vnormb+sqrt(g*B[0]))

def Fhatn(U):
    Uhat = U.Other(bnd=Ubnd)
    return (0.5*F(U)+0.5*F(Uhat))*n + Fmax(U,Uhat)/2*(U-Uhat)

DG formulation

def DGBilinearForm(fes,F,Fhatn,Ubnd):
    a = BilinearForm(fes, nonassemble=True)
    a += SymbolicBFI ( (- InnerProduct(F(U),grad(V).trans)).Compile())
    a += SymbolicBFI ( (InnerProduct(Fhatn(U),V)).Compile(), element_boundary=True)
    return a

a = DGBilinearForm(fes,F,Fhatn,Ubnd)

Simple fix to deal with shocks: artificial diffusion:

from DGdiffusion import AddArtificialDiffusion

artvisc = Parameter(1.0)
if order > 0:
    AddArtificialDiffusion(a,Ubnd,artvisc,compile=True)

Visualization of solution quantities

gfu = GridFunction(fes)
gfh, gfhu, gfhv = gfu
gfvu = gfhu/gfh
gfvv = gfhv/gfh
momentum = CoefficientFunction((gfhu,gfhv))
velocity = CoefficientFunction((gfvu,gfvv))
gfu.Set(U0)
Draw(momentum,mesh,"mom")
Draw(velocity,mesh,"vel")
Draw(gfh,mesh,"h")

Implicit Euler time stepping

def TimeLoop(a,gfu,dt,T,nsamplings=100):
    #gfu.Set(U0)
    res = gfu.vec.CreateVector()
    fes = a.space
    t = 0
    nsteps = int(ceil(T/dt))
    i = 0
    with TaskManager():
        while t <= T - 0.5*dt:
            a.Apply(gfu.vec, res)
            fes.SolveM(rho=CoefficientFunction(1),vec=res)
            gfu.vec.data -= dt * res
            t += dt
            if (i+1) % int(nsteps/nsamplings) == 0:
                Redraw()
            i+=1
            print("\rt=",t,end="")
    Redraw()
TimeLoop(a,gfu,dt=0.0004,T=5)

t= 4.9999999999995155745 0.14480000000000004 0.19120000000000137 0.2412000000000028 0.2908000000000042 0.3384000000000056 0.39320000000000715 0.5536000000000043 0.6071999999999984 0.6479999999999939 0.6991999999999883 0.7523999999999824 0.8007999999999771 0.9027999999999659 1.010399999999954 1.0315999999999517 1.22119999999993081.26919999999992551.4979999999999003 1.5783999999998914 1.7043999999998776 2.2131999999998215 2.3095999999998112.837199999999753 3.027199999999732 3.124399999999721 3.171999999999716 3.807599999999646 4.036799999999621 4.131199999999611 4.368399999999585

You may play around with this example.

• change the artificial diffusion parameter: How does it influence the solution
• change boundary conditions: left boundary -> fixed height and non-reflecting
• change the initial heights
• introduce a (circular) obstacle below the dam
• Generate output to create a video: To this end take a look at the final unit of this section.

%%HTML
<video width="600" height="400" controls>
    <source src="../../_static/shallow2D.mov">
</video>