added a 2d shallow water example

pydemo/shallowWater/main.py

@@ -6,7 +6,7 @@ from dune.femdg import createFemDGSolver
 
 from shallowwater import leVeque as problem
 
-Model, initial, x0,x1,N, endTime, name = problem()
+Model, initial, x0,x1,N, endTime, name = problem(1)
 
 parameter.append({"fem.verboserank": -1})
 parameter.append("parameter")
@@ -19,6 +19,7 @@ saveStep = 0.01
 saveTime = saveStep
 polOrder = 2
 limiter  = "default"
+subsamp  = 0
 
 space = create.space("dgonb", grid, order=polOrder, dimrange=dimR)
 u_h = space.interpolate(Model.toCons(initial), name='u_h')
@@ -27,7 +28,7 @@ operator = createFemDGSolver( Model, space, limiter=limiter )
 
 operator.applyLimiter( u_h );
 print("number of elements: ",grid.size(0),flush=True)
-vtk = grid.writeVTK(name, subsampling=2, write=False,
+vtk = grid.writeVTK(name, subsampling=subsamp, write=False,
            celldata=[u_h],
            pointdata={"freeSurface":eta},
            cellvector={"velocity":v}

pydemo/shallowWater/shallowwater.py

 from ufl import *
 from dune.ufl import Space
 
-def ShallowWater(topography):
+def ShallowWater(topography,g):
     dim = 2
     space = Space(2,3)
     x = SpatialCoordinate(space.cell())
-    g = 9.81
     class Model:
         dimension = dim+1
         def velo(U):
@@ -42,13 +41,24 @@ def ShallowWater(topography):
         boundary = {range(1,5): lambda t,x,u: u}
     return Model
 
-def leVeque():
-    cutoff = lambda x: conditional(2./5.<x[0],1,0)*conditional(x[0]<3./5.,1,0)
-    topography = lambda x: 1./4.*(cos(10*pi*(x[0]-0.5))+1)*cutoff(x)
-    space = Space(2,3)
-    x = SpatialCoordinate(space.cell())
-    initial = conditional(x[0]<0.1,1,conditional(x[0]<0.2,1.2,1))
-    return ShallowWater(topography),\
-           as_vector( [initial,0,0] ),\
-           [0, 0], [1, 0.25], [64, 16], 0.1,\
-           "leVeque"
+# example 5.1 and 7.1 from
+# https://www.sciencedirect.com/science/article/pii/S0021999198960582
+def leVeque(dim):
+    if dim == 1:
+        topography = lambda x: conditional(abs(x[0]-0.5)<0.1, 1./4.*(cos(10*pi*(x[0]-0.5))+1), 0)
+        space = Space(2,3)
+        x = SpatialCoordinate(space.cell())
+        initial = conditional(abs(x[0]-0.15)<0.05,1.2,1)
+        return ShallowWater(topography,1),\
+               as_vector( [initial,0,0] ),\
+               [0, 0], [1, 0.25], [64, 16], 0.7,\
+               "leVeque1D"
+    else:
+        topography = lambda x: 0.8*exp(-5*(x[0]-0.9)**2-50*(x[1]-0.5)**2)
+        space = Space(2,3)
+        x = SpatialCoordinate(space.cell())
+        initial = conditional(abs(x[0]-0.1)<0.05,1.01,1)
+        return ShallowWater(topography,1),\
+               as_vector( [initial,0,0] ),\
+               [0, 0], [2, 1], [80, 40], 1.8,\
+               "leVeque2D"