improve the testing environment a bit more and some further cleanup

Also make vortex euler problem work

pydemo/CMakeLists.txt

+dune_add_subdirs(scalar)
+dune_add_subdirs(shallowWater)
 dune_add_subdirs(euler)
-
-set(uflfiles
-  advection.ufl
-  diffusion.ufl
-  advectiondiffusion.ufl
-)
-
-dune_fem_add_elliptic_models(${uflfiles})
-
-set(headers)
-foreach(uflfile ${uflfiles})
-  get_filename_component(base ${uflfile} NAME_WE)
-  list(APPEND headers ${base}.hh)
-endforeach()
-
-add_custom_target(models DEPENDS ${headers})

pydemo/cpptest/CMakeLists.txt

+set(uflfiles
+  advection.ufl
+  diffusion.ufl
+  advectiondiffusion.ufl
+)
+
+dune_fem_add_elliptic_models(${uflfiles})
+
+set(headers)
+foreach(uflfile ${uflfiles})
+  get_filename_component(base ${uflfile} NAME_WE)
+  list(APPEND headers ${base}.hh)
+endforeach()
+
+add_custom_target(models DEPENDS ${headers})

pydemo/euler/__init__.py

+from .euler import problems

pydemo/euler/euler.py

@@ -63,7 +63,7 @@ def CompressibleEulerSlip(dim, gamma):
     return Model
 
 def riemanProblem(Model,x,x0,UL,UR):
-    return Model.toCons( conditional(x<x0,UL,UR) )
+    return Model.toCons( conditional(x<x0,as_vector(UL),as_vector(UR)) )
 
 # TODO Add exact solution where available (last argument)
 def constant(dim,gamma):
@@ -75,10 +75,8 @@ def sod(dim=2,gamma=1.4):
     space = Space(dim,dim+2)
     x = SpatialCoordinate(space.cell())
     Model = CompressibleEulerDirichlet(dim,gamma)
-    return  Model,\
-            riemanProblem( Model, x[0], 0.5,
-                          CompressibleEuler(dim,gamma).toCons([1,0,0,1]),
-                          CompressibleEuler(dim,gamma).toCons([0.125,0,0,0.1])),\
+    return Model,\
+           riemanProblem( Model, x[0], 0.5, [1,0,0,1], [0.125,0,0,0.1]),\
            [0, 0], [1, 0.25], [64, 16], 0.15,\
            "sod", None
 def radialSod1(dim=2,gamma=1.4):
@@ -87,8 +85,7 @@ def radialSod1(dim=2,gamma=1.4):
     Model = CompressibleEulerDirichlet(dim,gamma)
     return Model,\
            riemanProblem( Model, sqrt(dot(x,x)), 0.3,
-                          CompressibleEuler(dim,gamma).toCons([1,0,0,1]),
-                          CompressibleEuler(dim,gamma).toCons([0.125,0,0,0.1])),\
+                          [1,0,0,1], [0.125,0,0,0.1]),\
            [-0.5, -0.5], [0.5, 0.5], [20, 20], 0.25,\
            "radialSod1", None
 def radialSod1Large(dim=2,gamma=1.4):
@@ -97,8 +94,7 @@ def radialSod1Large(dim=2,gamma=1.4):
     Model = CompressibleEulerDirichlet(dim,gamma)
     return Model,\
            riemanProblem( Model, sqrt(dot(x,x)), 0.3,
-                          CompressibleEuler(dim,gamma).toCons([1,0,0,1]),
-                          CompressibleEuler(dim,gamma).toCons([0.125,0,0,0.1])),\
+                          [1,0,0,1], [0.125,0,0,0.1]),\
            [-1.5, -1.5], [1.5, 1.5], [60, 60], 0.5,\
            "radialSod1Large", None
 def radialSod2(dim=2,gamma=1.4):
@@ -107,8 +103,7 @@ def radialSod2(dim=2,gamma=1.4):
     Model = CompressibleEulerNeuman(dim,gamma)
     return Model,\
            riemanProblem( Model, sqrt(dot(x,x)), 0.3,
-                          CompressibleEuler(dim,gamma).toCons([0.125,0,0,0.1]),
-                          CompressibleEuler(dim,gamma).toCons([1,0,0,1])),\
+                          [0.125,0,0,0.1], [1,0,0,1]),\
            [-0.5, -0.5], [0.5, 0.5], [20, 20], 0.25,\
            "radialSod2", None
 def radialSod3(dim=2,gamma=1.4):
@@ -117,8 +112,7 @@ def radialSod3(dim=2,gamma=1.4):
     Model = CompressibleEulerSlip(dim,gamma)
     return Model,\
            riemanProblem( Model, sqrt(dot(x,x)), 0.3,
-                          CompressibleEuler(dim,gamma).toCons([1,0,0,1]),
-                          CompressibleEuler(dim,gamma).toCons([0.125,0,0,0.1])),\
+                          [1,0,0,1], [0.125,0,0,0.1]),\
            [-0.5, -0.5], [0.5, 0.5], [20, 20], 0.5,\
            "radialSod3", None
 
@@ -138,16 +132,16 @@ def vortex(dim=2,gamma=1.4):
     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)
+    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)
     Model = CompressibleEulerDirichlet(dim,gamma)
     return Model,\
            Model.toCons( as_vector( [rho,u,v,p] )),\
-           [-2, 2], [2, 2], [64, 64], 100,\
+           [-10, -10], [10, 10], [20, 20], 100,\
            "vortex", None
 
 problems = [sod, radialSod1, radialSod2, radialSod3,\

pydemo/euler/main.py

@@ -14,6 +14,7 @@ parameter.append("parameter")
 
 primitive=lambda Model,uh: {"pressure":Model.toPrim(uh)[2]}
 
-run(*problem(dim,gamma),
+run(*problem(),
         startLevel=0, polOrder=2, limiter="default",
-        primitive=primitive, saveStep=0.01, subsamp=2)
+        primitive=primitive, saveStep=0.1, subsamp=2,
+        dt=None)

pydemo/fulltest.py

+from dune.fem import parameter
+from dune.femdg.testing import run
+
+import scalar, shallowWater, euler
+
+parameter.append({"fem.verboserank": -1})
+parameter.append("parameter")
+
+for p in shallowWater.problems + euler.problems:
+    run(*p(), startLevel=0, polOrder=2, limiter="default",
+              primitive=None, saveStep=None, subsamp=0)

pydemo/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/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/parameter

+# LIMITER SETTINGS
+#-----------------
+# 0 = only dg solution | 1 = only reconstruction | 2 = both
+femdg.limiter.admissiblefunctions: 1
+# tolerance for shock indicator
+femdg.limiter.tolerance: 1 # 1e-10
+# 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/scalar/__init__.py

+from .scalar import problems

pydemo/scalar/scalar.py

@@ -83,6 +83,7 @@ class Burgers1D:
 
 space = Space(2,1)
 x = SpatialCoordinate(space.cell())
+
 def riemanProblem(x,x0,UL,UR):
     return conditional(x<x0,as_vector(UL),as_vector(UR))
 
@@ -141,3 +142,7 @@ def burgersVW():
 
 def burgersStationaryShock():
     return Burgers1D, riemanProblem(x[0],0,[1],[-1])
+
+problems = [ constantTransport, shockTransport, sinProblem,\
+             sinTransportProblem, pulse, diffusivePulse,\
+             burgersShock, burgersVW, burgersStationaryShock ]

pydemo/shallowWater/__init__.py

+from .shallowwater import problems

pydemo/shallowWater/shallowwater.py

@@ -61,3 +61,7 @@ def leVeque(dim):
                Model.toCons(as_vector( [initial,0,0] )),\
                [0, 0], [2, 1], [80, 40], 1.8,\
                "leVeque2D", None
+
+leVeque1D = lambda: leVeque(1)
+leVeque2D = lambda: leVeque(2)
+problems = [leVeque1D,leVeque2D]

python/dune/femdg/testing.py

@@ -8,7 +8,8 @@ from ufl import dot, SpatialCoordinate
 
 def run(Model, initial, x0,x1,N, endTime, name, exact,
         polOrder, limiter="default", startLevel=0,
-        primitive=None, saveStep=None, subsamp=0):
+        primitive=None, saveStep=None, subsamp=0,
+        dt=None):
     grid     = structuredGrid(x0,x1,N)
     grid.hierarchicalGrid.globalRefine(startLevel)
     dimR     = Model.dimension
@@ -19,7 +20,6 @@ def run(Model, initial, x0,x1,N, endTime, name, exact,
     space = create.space("dgonb", grid, order=polOrder, dimrange=dimR)
     u_h = space.interpolate(initial, name='u_h')
     operator = createFemDGSolver( Model, space, limiter=limiter )
-
     operator.applyLimiter( u_h );
     print("number of elements: ",grid.size(0),flush=True)
     if saveStep is not None:
@@ -38,18 +38,22 @@ def run(Model, initial, x0,x1,N, endTime, name, exact,
             velo[0].setConstant("time",[t])
         except:
             pass
-        vtk.write(name, count, outputType=OutputType.appendedraw)
+        vtk.write(name, count) # , outputType=OutputType.appendedraw)
     start = time.time()
     tcount = 0
     while t < endTime:
+        if dt is not None:
+            operator.setTimeStepSize(dt)
         operator.solve(target=u_h)
+        dt = operator.deltaT()
         if math.isnan( u_h.scalarProductDofs( u_h ) ):
+            grid.writeVTK(name, subsampling=subsamp, celldata=[u_h])
             print('ERROR: dofs invalid t =', t)
+            print('[',tcount,']','dt = ', dt, 'time = ',t, 'count = ',count, flush=True )
             exit(0)
-        dt = operator.deltaT()
         t += dt
         tcount += 1
-        if True: # tcount%100 == 0:
+        if tcount%100 == 0:
             print('[',tcount,']','dt = ', dt, 'time = ',t, 'count = ',count, flush=True )
         if saveStep is not None and t > saveTime:
             count += 1
@@ -68,6 +72,11 @@ def run(Model, initial, x0,x1,N, endTime, name, exact,
         except:
             pass
         vtk.write(name, count, outputType=OutputType.appendedraw)
+    else:
+        if exact is not None:
+            grid.writeVTK(name, subsampling=subsamp, celldata=[u_h, exact(t)])
+        else:
+            grid.writeVTK(name, subsampling=subsamp, celldata=[u_h])
 
     if exact is not None:
         error = integrate( grid, dot(u_h-exact(t),u_h-exact(t)), order=5 )