simplified the main script for the shallow water example

pydemo/shallowWater/main.py

@@ -3,7 +3,6 @@ from dune.grid import structuredGrid, cartesianDomain
 from dune.fem import parameter
 import dune.create as create
 from dune.femdg import createFemDGSolver
-# from ufl import *
 
 from shallowwater import leVeque as problem
 
@@ -12,54 +11,41 @@ Model, initial, x0,x1,N, endTime, name = problem()
 parameter.append({"fem.verboserank": -1})
 parameter.append("parameter")
 
-grid = structuredGrid(x0,x1,N)
-dimR     = grid.dimension + 1
+grid     = structuredGrid(x0,x1,N)
+dimR     = Model.dimension
 t        = 0
 count    = 0
 saveStep = 0.01
 saveTime = saveStep
-
-def initialize(space):
-    return space.interpolate(Model.toCons(initial), name='u_h')
-
-def useODESolver(polOrder=2, limiter='default'):
-    global count, t, saveTime
-    polOrder = polOrder
-    space = create.space("dgonb", grid, order=polOrder, dimrange=dimR)
-    u_h = initialize(space)
-    eta, v = Model.toPrim(u_h)
-    operator = createFemDGSolver( Model, space, limiter=limiter )
-
-    operator.applyLimiter( u_h );
-    print("number of elements: ",grid.size(0),flush=True)
-    grid.writeVTK(name,
-        celldata=[u_h],
-        pointdata={"freeSurface":eta},
-        cellvector={"velocity":v},
-        number=count, subsampling=2)
-    start = time.time()
-    tcount = 0
-    while t < endTime:
-        operator.solve(target=u_h)
-        dt = operator.deltaT()
-        t += dt
-        tcount += 1
-        if tcount%100 == 0:
-            print('[',tcount,']','dt = ', dt, 'time = ',t, 'count = ',count, flush=True )
-        if t > saveTime:
-            count += 1
-            grid.writeVTK(name,
-                celldata=[u_h],
-                pointdata={"freeSurface":eta},
-                cellvector={"velocity":v},
-                number=count, subsampling=2)
-            saveTime += saveStep
-    print("time loop:",time.time()-start)
-    print("number of time steps ", tcount)
-    grid.writeVTK(name,
-        celldata=[u_h],
-        pointdata={"freeSurface":eta},
-        cellvector={"velocity":v},
-        number=count, subsampling=2)
-
-useODESolver(2,'default')      # third order with limiter
+polOrder = 2
+limiter  = "default"
+
+space = create.space("dgonb", grid, order=polOrder, dimrange=dimR)
+u_h = space.interpolate(Model.toCons(initial), name='u_h')
+eta, v = Model.toPrim(u_h)
+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,
+           celldata=[u_h],
+           pointdata={"freeSurface":eta},
+           cellvector={"velocity":v}
+        )
+vtk.write(name, count)
+start = time.time()
+tcount = 0
+while t < endTime:
+    operator.solve(target=u_h)
+    dt = operator.deltaT()
+    t += dt
+    tcount += 1
+    if tcount%100 == 0:
+        print('[',tcount,']','dt = ', dt, 'time = ',t, 'count = ',count, flush=True )
+    if t > saveTime:
+        count += 1
+        vtk.write(name, count)
+        saveTime += saveStep
+print("time loop:",time.time()-start)
+print("number of time steps ", tcount)
+vtk.write(name, count)

pydemo/shallowWater/shallowwater.py

@@ -7,7 +7,7 @@ def ShallowWater(topography):
     x = SpatialCoordinate(space.cell())
     g = 9.81
     class Model:
-        dimension = dim+2
+        dimension = dim+1
         def velo(U):
             return as_vector( [U[i]/U[0] for i in range(1,dim+1)] )
         def toCons(U,x=x):
@@ -38,8 +38,8 @@ def ShallowWater(topography):
             return (hL - hR)/(0.5*(hL + hR))
         def S_ns(t,x,U,DU): # or S_s for a stiff source
             return as_vector( [0, *(-U[0]*g*grad(topography(x))) ])
-        boundary = {range(1,5): lambda t,x,u,n: Model.F_c(t,x,u)*n}
-        # boundary = {range(1,5): lambda t,x,u: u}
+        # boundary = {range(1,5): lambda t,x,u,n: Model.F_c(t,x,u)*n}
+        boundary = {range(1,5): lambda t,x,u: u}
     return Model
 
 def leVeque():
@@ -48,7 +48,7 @@ def leVeque():
     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] ) ,\
+    return ShallowWater(topography),\
+           as_vector( [initial,0,0] ),\
            [0, 0], [1, 0.25], [64, 16], 0.1,\
            "leVeque"