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

3.4 A Nonlinear conservation law: 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) \end{align*}

\begin{align*} \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) \end{align*}

\begin{align*} \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} \}\]
[1]:

from ngsolve import * # sloppy
import netgen.gui

The dam break problem (geometry)

alt

[2]:

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]

[3]:

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

[4]:

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

 Generate Mesh from spline geometry
 Boundary mesh done, np = 59
 CalcLocalH: 59 Points 0 Elements 0 Surface Elements
 Meshing domain 1 / 2
 load internal triangle rules
 Surface meshing done
 Meshing domain 2 / 2
 Surface meshing done
 Edgeswapping, topological
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Update mesh topology
 Update clusters
 Curve elements, order = 3
 Update clusters
 Curving edges
 Curving faces

Vectorial (dim=3) approximation space:

[5]:

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

initial and boundary conditions

[6]:

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}\]
[7]:

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})\]
[8]:

n = specialcf.normal(mesh.dim)

def Max(u,v):
    return IfPos(u-v,u,v)
def Fmax(A,B): # max. characteristic speed:
    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): # numerical flux
    Uhat = U.Other(bnd=Ubnd)
    return (0.5*F(U)+0.5*F(Uhat))*n + Fmax(U,Uhat)/2*(U-Uhat)

DG formulation

We recall that a BilinearForm is allowed to be nonlinear in the first argument.

[9]:

def DGBilinearForm(fes,F,Fhatn,Ubnd):
    a = BilinearForm(fes, nonassemble=True)
    a += - InnerProduct(F(U),Grad(V)) * dx
    a += InnerProduct(Fhatn(U),V) * dx(element_boundary=True)
    return a

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

Simple fix to deal with shocks: artificial diffusion:

[10]:

from DGdiffusion import AddArtificialDiffusion

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

Compiled CF:
N5ngfem28ParameterCoefficientFunctionIdEE
N5ngfem24ScaleCoefficientFunctionE
Z25ExportCoefficientFunctionRN8pybind117module_EE10MeshSizeCF
N5ngfem13cl_BinaryOpCFINS_11GenericMultEEE
N5ngfem27ConstantCoefficientFunctionE
N5ngfem13cl_BinaryOpCFINS_10GenericDivEEE
N5ngfem13ProxyFunctionE
N5ngfem13ProxyFunctionE
N5ngfem31T_MultVecVecCoefficientFunctionILi6EEE
N5ngfem13cl_BinaryOpCFINS_11GenericMultEEE
inputs =
0:
1: 0
2:
3: 1 2
4:
5: 3 4
6:
7:
8: 6 7
9: 5 8

Compiled CF:
N5ngfem28ParameterCoefficientFunctionIdEE
N5ngfem24ScaleCoefficientFunctionE
Z25ExportCoefficientFunctionRN8pybind117module_EE10MeshSizeCF
N5ngfem13cl_BinaryOpCFINS_11GenericMultEEE
N5ngfem27ConstantCoefficientFunctionE
N5ngfem13cl_BinaryOpCFINS_10GenericDivEEE
N5ngfem13ProxyFunctionE
N5ngfem28TransposeCoefficientFunctionE
N5ngfem24ScaleCoefficientFunctionE
N5ngfem13ProxyFunctionE
N5ngfem28TransposeCoefficientFunctionE
N5ngfem24ScaleCoefficientFunctionE
N5ngfem13cl_BinaryOpCFINS_11GenericPlusEEE
N5ngfem17cl_NormalVectorCFILi2EEE
N5ngfem29MultMatVecCoefficientFunctionE
N5ngfem13ProxyFunctionE
N5ngfem13ProxyFunctionE
N5ngfem13cl_BinaryOpCFINS_12GenericMinusEEE
N5ngfem31T_MultVecVecCoefficientFunctionILi3EEE
N5ngfem24ScaleCoefficientFunctionE
N5ngfem13ProxyFunctionE
N5ngfem28TransposeCoefficientFunctionE
N5ngfem24ScaleCoefficientFunctionE
N5ngfem13ProxyFunctionE
N5ngfem28TransposeCoefficientFunctionE
N5ngfem24ScaleCoefficientFunctionE
N5ngfem13cl_BinaryOpCFINS_11GenericPlusEEE
N5ngfem29MultMatVecCoefficientFunctionE
N5ngfem13ProxyFunctionE
N5ngfem13ProxyFunctionE
N5ngfem13cl_BinaryOpCFINS_12GenericMinusEEE
N5ngfem31T_MultVecVecCoefficientFunctionILi3EEE
N5ngfem13cl_BinaryOpCFINS_12GenericMinusEEE
N5ngfem27ConstantCoefficientFunctionE
N5ngfem13cl_BinaryOpCFINS_10GenericDivEEE
N5ngfem13cl_BinaryOpCFINS_12GenericMinusEEE
N5ngfem13ProxyFunctionE
N5ngfem13cl_BinaryOpCFINS_12GenericMinusEEE
N5ngfem31T_MultVecVecCoefficientFunctionILi3EEE
N5ngfem13cl_BinaryOpCFINS_11GenericMultEEE
N5ngfem13cl_BinaryOpCFINS_11GenericPlusEEE
N5ngfem13cl_BinaryOpCFINS_11GenericMultEEE
inputs =
0:
1: 0
2:
3: 1 2
4:
5: 3 4
6:
7: 6
8: 7
9:
10: 9
11: 10
12: 8 11
13:
14: 12 13
15:
16:
17: 15 16
18: 14 17
19: 18
20:
21: 20
22: 21
23:
24: 23
25: 24
26: 22 25
27: 26 13
28:
29:
30: 28 29
31: 27 30
32: 19 31
33:
34: 33 2
35: 28 29
36:
37: 15 36
38: 35 37
39: 34 38
40: 32 39
41: 5 40

Visualization of solution quantities

[11]:

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

setvalues element 206/206

Explicit Euler time stepping

[12]:

def TimeLoop(a,gfu,dt,T,nsamplings=100):
    #gfu.Set(U0)
    res = gfu.vec.CreateVector()
    fes = a.space
    t = 0; i = 0
    nsteps = int(ceil(T/dt))
    invma = fes.Mass(1).Inverse() @ a.mat
    with TaskManager():
        while t <= T - 0.5*dt:
            res = invma * gfu.vec
            gfu.vec.data -= dt * res
            t += dt
            if (i+1) % int(nsteps/nsamplings) == 0:
                Redraw()
            i+=1
            print("\rt = {:.10}".format(t),end="")
    Redraw()

[13]:

TimeLoop(a,gfu,dt=0.0004,T=3)

t = 3.0

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.

[14]:

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