first overhaul of python script structure - adding a testing script to

python/dune/femdg
that can be used for testing the different models and problems posed on
Cartesian meshes.

pydemo/euler/euler.py

@@ -3,6 +3,7 @@ from dune.ufl import Space
 
 def CompressibleEuler(dim, gamma):
     class Model:
+        dimension = dim+2
         def velo(U):
             return as_vector( [U[i]/U[0] for i in range(1,dim+1)] )
         def rhoeps(U):
@@ -64,11 +65,12 @@ def CompressibleEulerSlip(dim, gamma):
 def riemanProblem(x,x0,UL,UR):
     return conditional(x<x0,UL,UR)
 
+# TODO Add exact solution where available (last argument)
 def constant(dim,gamma):
     return CompressibleEulerDirichlet(dim,gamma) ,\
            as_vector( [0.1,0.,0.,0.1] ),\
            [-1, 0], [1, 0.1], [50, 5], 0.1,\
-           "constant"
+           "constant", None
 def sod(dim,gamma):
     space = Space(dim,dim+2)
     x = SpatialCoordinate(space.cell())
@@ -77,7 +79,7 @@ def sod(dim,gamma):
                           CompressibleEuler(dim,gamma).toCons([1,0,0,1]),
                           CompressibleEuler(dim,gamma).toCons([0.125,0,0,0.1])),\
            [0, 0], [1, 0.25], [64, 16], 0.15,\
-           "sod"
+           "sod", None
 def radialSod1(dim,gamma):
     space = Space(dim,dim+2)
     x = SpatialCoordinate(space.cell())
@@ -86,7 +88,7 @@ def radialSod1(dim,gamma):
                           CompressibleEuler(dim,gamma).toCons([1,0,0,1]),
                           CompressibleEuler(dim,gamma).toCons([0.125,0,0,0.1])),\
            [-0.5, -0.5], [0.5, 0.5], [20, 20], 0.25,\
-           "radialSod1"
+           "radialSod1", None
 def radialSod1Large(dim,gamma):
     space = Space(dim,dim+2)
     x = SpatialCoordinate(space.cell())
@@ -95,7 +97,7 @@ def radialSod1Large(dim,gamma):
                           CompressibleEuler(dim,gamma).toCons([1,0,0,1]),
                           CompressibleEuler(dim,gamma).toCons([0.125,0,0,0.1])),\
            [-1.5, -1.5], [1.5, 1.5], [60, 60], 0.5,\
-           "radialSod1Large"
+           "radialSod1Large", None
 def radialSod2(dim,gamma):
     space = Space(dim,dim+2)
     x = SpatialCoordinate(space.cell())
@@ -104,7 +106,7 @@ def radialSod2(dim,gamma):
                           CompressibleEuler(dim,gamma).toCons([0.125,0,0,0.1]),
                           CompressibleEuler(dim,gamma).toCons([1,0,0,1])),\
            [-0.5, -0.5], [0.5, 0.5], [20, 20], 0.25,\
-           "radialSod2"
+           "radialSod2", None
 def radialSod3(dim,gamma):
     space = Space(dim,dim+2)
     x = SpatialCoordinate(space.cell())
@@ -113,7 +115,7 @@ def radialSod3(dim,gamma):
                           CompressibleEuler(dim,gamma).toCons([1,0,0,1]),
                           CompressibleEuler(dim,gamma).toCons([0.125,0,0,0.1])),\
            [-0.5, -0.5], [0.5, 0.5], [20, 20], 0.5,\
-           "radialSod3"
+           "radialSod3", None
 
 def leVeque(dim,gamma):
     space = Space(dim,dim+2)
@@ -122,4 +124,21 @@ def leVeque(dim,gamma):
     return CompressibleEulerDirichlet(dim,gamma),\
            as_vector( [initial,0,0,initial] ),\
            [0, 0], [1, 0.25], [64, 16], 0.7,\
-           "leVeque1D"
+           "leVeque1D", None
+
+def vortex(dim,gamma):
+    S = 13.5    # strength of vortex
+    R = 1.5     # radius of vortex
+    M = 0.4     # Machnumber
+
+    space = Space(dim,dim+2)
+    x = SpatialCoordinate(space.cell())
+    f = (1 - x[0]*x[0] - x[1]*x[1])/(2*R*R)
+    rho = pow(1 - S*S*M*M*(gamma - 1)*exp(2*f)/(8*pi*pi), 1/(gamma - 1))
+    u =     S*x[1]*exp(f)/(2*pi*R)
+    v = 1 - S*x[0]*exp(f)/(2*pi*R)
+    p = rho / (gamma*M*M)
+    return CompressibleEulerDirichlet(dim,gamma),\
+           as_vector( [rho,u,v,p] ),\
+           [-2, 2], [2, 2], [64, 64], 100,\
+           "vortex", None

pydemo/euler/main.py

+from dune.fem import parameter
+from dune.femdg.testing import run
+
+from euler import sod as problem
+# from euler import vortex as problem
+# from euler import leVeque as problem
+# from euler import radialSod3 as problem
+
+dim = 2
+gamma = 1.4
+
+parameter.append({"fem.verboserank": -1})
+parameter.append("parameter")
+
+primitive=lambda Model,uh: {"freesurface":Model.toPrim(uh)[0]}
+
+run(*problem(dim,gamma),
+        startLevel=0, polOrder=2, limiter="default",
+        primitive=primitive, saveStep=0.01, subsamp=2)

pydemo/euler/shock_tube.py

-import time
-from dune.grid import structuredGrid, cartesianDomain
-from dune.fem import parameter
-import dune.create as create
-from dune.models.elliptic.formfiles import loadModels
-from llf import NumFlux
-from dune.femdg import createFemDGSolver
-from ufl import *
-
-gamma = 1.4
-dim = 2
-
-from euler import sod as problem
-# from euler import leVeque as problem
-# from euler import radialSod3 as problem
-
-Model, initial, x0,x1,N, endTime, name = problem(dim,gamma)
-
-parameter.append({"fem.verboserank": -1})
-parameter.append("parameter")
-
-grid = structuredGrid(x0,x1,N)
-# grid = create.grid("ALUSimplex", cartesianDomain(x0,x1,N))
-dimR     = grid.dimension + 2
-t        = 0
-dt       = 1e-3
-count    = 0
-saveStep = 0.01
-saveTime = saveStep
-
-def initialize(space):
-    return space.interpolate(initial, name='u_h')
-    if space.order == 0:
-        return space.interpolate(initial, name='u_h')
-    else:
-        lagOrder = 1 # space.order
-        spacelag = create.space("lagrange", space.grid, order=lagOrder, dimrange=space.dimRange)
-        u_h = spacelag.interpolate(initial, name='tmp')
-        return space.interpolate(u_h, name='u_h')
-
-def useGalerkinOp():
-    global count, t, dt, saveTime
-    # a full fv implementation using UFL and the galerkin operator
-    # some sign error or something else leading to wrong results
-    # still needs fixing
-    spaceName = "finitevolume"
-    space = create.space(spaceName, grid, dimrange=dimR)
-
-    n = FacetNormal(space.cell())
-
-    u_h   = space.interpolate(initial, name='u_h')
-    u_h_n = u_h.copy(name="previous")
-
-    fullModel = inner( Model.F_c(u), grad(v) ) * dx -\
-                inner( NumFlux(Model,u,n,dt), v('+')-v('-')) * dS -\
-                inner( Model.F_c(u)*n, v) * ds
-    operator = create.operator("galerkin",
-                 inner(u-u_h_n,v)*dx == dt*replace(fullModel,{u:u_h_n}),
-                 space )
-
-    operator.model.dt = 1e-5
-
-    grid.writeVTK(name, pointdata=[u_h], number=count)
-    start = time.time()
-    while t < endTime:
-        u_h_n.assign(u_h)
-        operator.solve(target=u_h)
-        t += operator.model.dt
-        if t > saveTime:
-            count += 1
-            grid.writeVTK(name, pointdata=[u_h], number=count)
-            saveTime += saveStep
-    print("time loop:",time.time()-start)
-    grid.writeVTK(name, pointdata=[u_h], number=count)
-
-def useODESolver(polOrder=2, limiter='default'):
-    global count, t, dt, saveTime
-    polOrder = polOrder
-    if False:
-        # needs change in dune/fem-dg/operator/dg/passtraits.hh
-        space = create.space("dglegendre", grid, order=polOrder, dimrange=dimR, hierarchical=False)
-    else:
-        space = create.space("dgonb", grid, order=polOrder, dimrange=dimR)
-    u_h = initialize(space)
-    rho, v, p = Model.toPrim(u_h)
-    operator = createFemDGSolver( Model, space, limiter=limiter )
-    # operator.setTimeStepSize(dt)
-
-    operator.applyLimiter( u_h );
-    print("number of elements: ",grid.size(0),flush=True)
-    grid.writeVTK(name,
-        pointdata=[u_h],
-        # celldata={"density":rho, "pressure":p}, # bug: density not shown correctly
-        celldata={"pressure":p, "maxLambda":Model.maxLambda(0,0,u_h,as_vector([1,0]))},
-        cellvector={"velocity":v},
-        number=count, subsampling=2)
-    start = time.time()
-    tcount = 0
-    while t < endTime:
-        operator.solve(target=u_h)
-        # operator.applyLimiter( 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,
-                pointdata=[u_h],
-                celldata={"pressure":p, "maxLambda":Model.maxLambda(0,0,u_h,as_vector([1,0]))},
-                cellvector={"velocity":v},
-                number=count, subsampling=2)
-            saveTime += saveStep
-    print("time loop:",time.time()-start)
-    print("number of time steps ", tcount)
-    grid.writeVTK(name,
-        pointdata=[u_h],
-        celldata={"pressure":p, "maxLambda":Model.maxLambda(0,0,u_h,as_vector([1,0]))},
-        cellvector={"velocity":v},
-        number=count, subsampling=2)
-
-if True:
-    # grid = structuredGrid(x0,x1,N)
-    grid = create.grid("ALUCube", cartesianDomain(x0,x1,N))
-    grid.hierarchicalGrid.globalRefine(1)
-    # grid = create.view("adaptive", grid)
-    useODESolver(2,'default')      # third order with limiter
-elif False:
-    N = [n*10 for n in N]
-    grid = structuredGrid(x0,x1,N)
-    # grid = create.grid("ALUSimplex", cartesianDomain(x0,x1,N))
-    useODESolver(0,None)           # FV scheme
-elif False:
-    N = [n*10 for n in N]
-    grid = structuredGrid(x0,x1,N)
-    # grid = create.grid("ALUSimplex", cartesianDomain(x0,x1,N))
-    useODESolver(0,'default')      # FV scheme with limiter

pydemo/scalar/main.py

+from dune.fem import parameter
+from dune.femdg.testing import run
+
+# from scalar import shockTransport as problem
+# from scalar import sinProblem as problem
+from scalar import sinTransportProblem as problem
+# from scalar import pulse as problem
+# from scalar import diffusivePulse as problem
+
+parameter.append({"fem.verboserank": 0})
+parameter.append("parameter")
+
+run(*problem(),
+        startLevel=0, polOrder=2, limiter=None,
+        primitive=None, saveStep=0.1, subsamp=0)

pydemo/scalar/scalar.py

 from ufl import *
 from dune.ufl import Space
 class Transport1D:
+    dimension = 1
     name = 'transport'
 
+    def velo(U):
+        return as_vector( [1,0] )
+    def toCons(U):
+        return U
     def F_c(t,x,U):
         return as_matrix( [ [U[0], 0] ] )
 
@@ -11,7 +16,7 @@ class Transport1D:
     def maxLambda(t,x,U,n):
         return abs(n[0])
     def velocity(t,x,U):
-        return as_vector( [1,0] )
+        return Transport1D.velo(U)
     def physical(U):
         return 1
     def jump(U,V):
@@ -20,6 +25,9 @@ class Transport1D:
 def LinearAdvectionDiffusion1D(v,eps):
     if v is not None: v = as_vector(v)
     class Model:
+        dimension = 1
+        def toCons(U):
+            return U
         if v is not None and eps is not None:
             name = 'linearAD'
         elif v is not None:
@@ -27,13 +35,17 @@ def LinearAdvectionDiffusion1D(v,eps):
         else:
             name = 'linearD'
         if v is not None:
+            def velo(U):
+                return v
             def F_c(t,x,U):
                 return as_matrix( [[ *(v*U[0]) ]] )
             def maxLambda(t,x,U,n):
                 return abs(dot(v,n))
             def velocity(t,x,U):
-                return v
+                return Model.velo(U)
         if eps is not None:
+            def velo(U):
+                return as_vector([0,0])
             def F_v(t,x,U,DU):
                 return eps*DU
             # TODO: this method is used in an IMEX method allough the
@@ -68,8 +80,13 @@ def LinearAdvectionDiffusion1DDirichlet(v,eps,bnd):
 
 # burgers problems still need to be adapted to new API
 class Burgers1D:
+    dimension = 1
     name = 'burgers'
+    def toCons(U):
+        return U
 
+    def velo(U):
+        return as_vector( [U[0],0] )
     def F_c(t,x,U):
         return as_matrix( [ [U[0]*U[0]/2, 0] ] )
 
@@ -81,7 +98,7 @@ class Burgers1D:
     def maxLambda(t,x,U,n):
         return abs(U[0]*n[0])
     def velocity(t,x,U):
-        return as_vector( [U[0],0] )
+        return Burgers1D.velo(U)
     def physical(U):
         return 1
     def jump(U,V):

pydemo/scalar/shock_tube.py

-import time, math
-from dune.grid import structuredGrid, cartesianDomain
-from dune.fem import parameter
-from dune.fem.function import integrate
-import dune.create as create
-from dune.models.elliptic.formfiles import loadModels
-from dune.femdg import createFemDGSolver
-from ufl import dot
-
-# from scalar import shockTransport as problem
-# from scalar import sinProblem as problem
-from scalar import sinTransportProblem as problem
-# from scalar import pulse as problem
-# from scalar import diffusivePulse as problem
-
-Model, initial, x0,x1,N, endTime, name, exact = problem()
-
-parameter.append({"fem.verboserank": 0})
-parameter.append("parameter")
-
-grid = structuredGrid(x0,x1,N)
-
-# grid = create.grid("ALUSimplex", cartesianDomain(x0,x1,N))
-dimR     = 1
-t        = 0
-dt       = 0.001
-saveStep = 0.1
-saveTime = saveStep
-count    = 0
-
-def useODESolver():
-    global count, t, dt, saveTime
-
-    spaceName = "dgonb"
-    polOrder = 2
-    space = create.space(spaceName, grid, order=polOrder, dimrange=dimR)
-    u_h   = space.interpolate(initial, name='u_h')
-    operator = createFemDGSolver( Model, space, limiter=None, diffusionScheme="cdg2" )
-
-    #operator.setTimeStepSize(dt)
-
-    start = time.time()
-    #grid.writeVTK(name, celldata=[u_h], number=count, subsampling=2)
-    grid.writeVTK(name, celldata=[u_h], number=count)
-    print('dt = ', dt, 'time = ',t, 'count = ',count, flush=True )
-    while t < endTime:
-        operator.solve(target=u_h)
-        if math.isnan( u_h.scalarProductDofs( u_h ) ):
-            print('ERROR: dofs invalid t =', t)
-            exit(0)
-        dt = operator.deltaT()
-        t += dt
-        print('dt = ', dt, 'time = ',t, 'count = ',count, flush=True )
-        if t > saveTime:
-            count += 1
-            grid.writeVTK(name, celldata=[u_h], number=count)
-            saveTime += saveStep
-    grid.writeVTK(name, celldata=[u_h], number=count)
-    print("time loop:",time.time()-start)
-    if exact is not None:
-        error = integrate( grid, dot(u_h-exact(t),u_h-exact(t)), order=5 )
-        print("error:", math.sqrt(error) )
-
-useODESolver()

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 dune.femdg.testing import run
 
 from shallowwater import leVeque as problem
 
-Model, initial, x0,x1,N, endTime, name = problem(1)
-
 parameter.append({"fem.verboserank": -1})
 parameter.append("parameter")
 
-grid     = structuredGrid(x0,x1,N)
-# grid.hierarchicalGrid.globalRefine(3)
-dimR     = Model.dimension
-t        = 0
-count    = 0
-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')
-eta, v = Model.toPrim(u_h)
-operator = createFemDGSolver( Model, space, limiter=limiter )
+primitive=lambda Model,uh: {"freesurface":Model.toPrim(uh)[0]}
 
-operator.applyLimiter( u_h );
-print("number of elements: ",grid.size(0),flush=True)
-vtk = grid.writeVTK(name, subsampling=subsamp, 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)
+run(*problem(1),
+        startLevel=0, polOrder=2, limiter="default",
+        primitive=primitve, saveStep=0.1, subsamp=0)

pydemo/shallowWater/shallowwater.py

@@ -33,8 +33,8 @@ def ShallowWater(topography,g):
             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
 
 # example 5.1 and 7.1 from
@@ -48,7 +48,7 @@ def leVeque(dim):
         return ShallowWater(topography,1),\
                as_vector( [initial,0,0] ),\
                [0, 0], [1, 0.25], [64, 16], 0.7,\
-               "leVeque1D"
+               "leVeque1D", None
     else:
         topography = lambda x: 0.8*exp(-5*(x[0]-0.9)**2-50*(x[1]-0.5)**2)
         space = Space(2,3)
@@ -57,4 +57,4 @@ def leVeque(dim):
         return ShallowWater(topography,1),\
                as_vector( [initial,0,0] ),\
                [0, 0], [2, 1], [80, 40], 1.8,\
-               "leVeque2D"
+               "leVeque2D", None

python/dune/femdg/CMakeLists.txt

 add_python_targets(fem-dg
   __init__
   _operators
+  testing
 )

python/dune/femdg/testing.py

+import time, math
+from dune.grid import structuredGrid
+import dune.create as create
+from dune.femdg import createFemDGSolver
+
+def run(Model, initial, x0,x1,N, endTime, name, exact,
+        polOrder, limiter="default", startLevel=0,
+        primitive, saveStep, subsamp=0):
+    grid     = structuredGrid(x0,x1,N)
+    grid.hierarchicGrid.globalRefine(startLevel)
+    dimR     = Model.dimension
+    t        = 0
+    count    = 0
+    saveTime = saveStep
+
+    space = create.space("dgonb", grid, order=polOrder, dimrange=dimR)
+    u_h = space.interpolate(Model.toCons(initial), name='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=subsamp, write=False,
+               celldata=[u_h],
+               pointdata=primitive(Model,u_h) if primitive else None,
+               cellvector={"velocity":Model.velo(u_h)}
+            )
+    vtk.write(name, count)
+    start = time.time()
+    tcount = 0
+    while t < endTime:
+        operator.solve(target=u_h)
+        if math.isnan( u_h.scalarProductDofs( u_h ) ):
+            print('ERROR: dofs invalid t =', t)
+            exit(0)
+        dt = operator.deltaT()
+        t += dt
+        tcount += 1
+        if True: # 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)
+
+    if exact is not None:
+        error = integrate( grid, dot(u_h-exact(t),u_h-exact(t)), order=5 )
+        print("error:", math.sqrt(error) )