added an example with a source term (shallow water)

Some issue with the boundary conditions need to be looked into

pydemo/shallowWater/main.py

+import time
+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
+
+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
+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

pydemo/shallowWater/paramBase

+# OMP THREADS 
+#------------
+fem.verboserank: 0
+# number of threads used in OMP program
+fem.parallel.numberofthreads: 1
+# write diagnostics file (
+# 0 don't, 1 only speedup file, 2 write all runfiles 
+# 3 only write 0, others at end, 4 all files at end for scaling) 
+fem.parallel.diagnostics: 1
+# if true non-blocking communication is enabled
+femdg.nonblockingcomm: true
+fem.threads.verbose: true
+# if true thread 0 only communicates data but does not computation
+fem.threads.communicationthread: false
+fem.threads.partitioningmethod: sfc

pydemo/shallowWater/paramSolver

+# SOLVER CONFIGURATION
+#---------------------
+
+fem.ode.odesolver: EX # ode solvers: EX, IM, IMEX
+# fem.ode.order: 3
+fem.ode.verbose: cfl # ode output: none, cfl, full
+fem.ode.cflincrease: 1.25
+fem.ode.miniterations: 95
+fem.ode.maxiterations: 105
+fem.ode.cflStart: 1.
+#fem.ode.cflMax: 5
+fem.timeprovider.factor: 0.45
+fem.timeprovider.updatestep: 1
+femdg.stepper.maxtimestep: 1e-3
+
+# parameter for the implicit solvers
+fem.solver.verbose: false
+fem.solver.gmres.restart: 15
+fem.solver.newton.verbose: false
+fem.solver.newton.linear.verbose: false
+fem.solver.newton.maxlineariterations: 1000
+fem.solver.newton.tolerance: 1e-10
+
+dgdiffusionflux.method: CDG2 # diffusion flux: CDG2, CDG, BR2, IP, NIPG, BO
+dgdiffusionflux.theoryparameters: 1 # scaling with theory parameters
+dgdiffusionflux.penalty: 0.
+dgdiffusionflux.liftfactor: 1.0

pydemo/shallowWater/parameter

+# LIMITER SETTINGS
+#-----------------
+# 0 = only dg solution | 1 = only reconstruction | 2 = both
+femdg.limiter.admissiblefunctions: 1
+# tolerance for shock indicator
+femdg.limiter.tolerance: 1
+# threshold for avoiding over-excessive limitation
+femdg.limiter.limiteps: 1e-8
+
+# GENERAL
+#--------
+paramfile: paramBase
+
+
+# choises are: LLF, HLL, HLLC, LLF2
+dgadvectionflux.method: LLF
+
+# SOLVER CONFIGURATION
+#---------------------
+fem.ode.odesolver: EX
+paramfile: paramSolver

pydemo/shallowWater/shallowwater.py

+from ufl import *
+from dune.ufl import Space
+
+def ShallowWater(topography):
+    dim = 2
+    space = Space(2,3)
+    x = SpatialCoordinate(space.cell())
+    g = 9.81
+    class Model:
+        dimension = dim+2
+        def velo(U):
+            return as_vector( [U[i]/U[0] for i in range(1,dim+1)] )
+        def toCons(U,x=x):
+            v = as_vector( [U[i] for i in range(1,dim+1)] )
+            return as_vector( [U[0]-topography(x), *(U[0]*v)] )
+        def toPrim(U,x=x):
+            return U[0]+topography(x), Model.velo(U)
+        def F_c(t,x,U):
+            assert dim==2
+            h = U[0]
+            v = Model.velo(U)
+            p = 0.5*g*h*h
+            return as_matrix([
+                    [h*v[0], h*v[1]],
+                    [h*v[0]*v[0] + p, h*v[0]*v[1]],
+                    [h*v[0]*v[1], h*v[1]*v[1] + p] ] )
+        def maxLambda(t,x,U,n):
+            h = U[0]
+            v = Model.velo(U)
+            return abs(dot(v,n)) + sqrt(g*h)
+        def velocity(t,x,U):
+            return Model.velo(U)
+        def physical(U):
+            return conditional( (U[0]>1e-8), 1, 0 )
+        def jump(U,V):
+            hL = U[0]
+            hR = V[0]
+            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}
+    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"