improve API for boundary conditions - use a single dictionary and allow to

provide tuples/lists or ranges of ids

pydemo/euler/euler.py

@@ -47,25 +47,18 @@ def CompressibleEuler(dim, gamma):
 
 def CompressibleEulerNeuman(dim, gamma):
     class Model(CompressibleEuler(dim,gamma)):
-        def outflowFlux(t,x,u,n):
-            return Model.F_c(t,x,u)*n
-        boundaryFlux = {}
-        for i in range(1,5): boundaryFlux.update( {i: outflowFlux} )
+        boundary = {range(1,5): lambda t,x,u,n: Model.F_c(t,x,u)*n}
     return Model
 def CompressibleEulerDirichlet(dim, gamma):
     class Model(CompressibleEuler(dim,gamma)):
-        def outflowValue(t,x,u):
-            return u
-        boundaryValue = {}
-        for i in range(1,5): boundaryValue.update( {i: outflowValue} )
+        boundary = {range(1,5): lambda t,x,u: u}
     return Model
 def CompressibleEulerSlip(dim, gamma):
     class Model(CompressibleEuler(dim,gamma)):
         def outflowFlux(t,x,u,n):
             _,_, p = CompressibleEuler(dim,gamma).toPrim(u)
             return as_vector([ 0, *(p*n), 0 ])
-        boundaryFlux = {}
-        for i in range(1,5): boundaryFlux.update( {i: outflowFlux} )
+        boundary = {range(1,5): outflowFlux}
     return Model
 
 def riemanProblem(x,x0,UL,UR):

pydemo/euler/paramBase

@@ -2,7 +2,7 @@
 #------------
 fem.verboserank: 0
 # number of threads used in OMP program
-fem.parallel.numberofthreads: 4
+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) 

pydemo/euler/paramSolver

@@ -2,7 +2,7 @@
 #---------------------
 
 fem.ode.odesolver: EX # ode solvers: EX, IM, IMEX
-fem.ode.order: 3
+# fem.ode.order: 3
 fem.ode.verbose: cfl # ode output: none, cfl, full
 fem.ode.cflincrease: 1.25
 fem.ode.miniterations: 95

pydemo/euler/shock_tube.py

@@ -9,8 +9,9 @@ from ufl import *
 
 gamma = 1.4
 dim = 2
-from euler import sod as problem
-# from euler import radialSod3 as problem
+
+# from euler import sod as problem
+from euler import radialSod3 as problem
 
 Model, initial, x0,x1,N, endTime, name = problem(dim,gamma)
 
@@ -27,6 +28,7 @@ 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:
@@ -83,7 +85,7 @@ def useODESolver(polOrder=2, limiter='default'):
     operator = createFemDGSolver( Model, space, limiter=limiter )
     # operator.setTimeStepSize(dt)
 
-    # operator.applyLimiter( u_h );
+    operator.applyLimiter( u_h );
     print("number of elements: ",grid.size(0),flush=True)
     grid.writeVTK(name,
         pointdata=[u_h],
@@ -101,9 +103,8 @@ def useODESolver(polOrder=2, limiter='default'):
         tcount += 1
         if tcount%100 == 0:
             print('[',tcount,']','dt = ', dt, 'time = ',t, 'count = ',count, flush=True )
-        if False: # t > saveTime:
+        if t > saveTime:
             count += 1
-            # rho, v, p = Model.toPrim(u_h) # is this needed - works for me # without and it should...
             grid.writeVTK(name,
                 pointdata=[u_h],
                 celldata={"pressure":p, "maxLambda":Model.maxLambda(0,0,u_h,as_vector([1,0]))},
@@ -121,6 +122,8 @@ def useODESolver(polOrder=2, limiter='default'):
 if True:
     # grid = structuredGrid(x0,x1,N)
     grid = create.grid("ALUCube", cartesianDomain(x0,x1,N))
+    # grid.hierarchicalGrid.globalRefine(2)
+    # grid = create.view("adaptive", grid)
     useODESolver(2,'default')      # third order with limiter
 elif False:
     N = [n*10 for n in N]

pydemo/scalar/paramSolver

@@ -15,16 +15,17 @@ fem.timeprovider.updatestep: 1
 # parameter for the implicit solvers
 # fem.differenceoperator.eps: 1e-12
 fem.solver.gmres.restart: 50
-fem.solver.newton.verbose: false
+fem.solver.newton.verbose: true
 fem.solver.newton.linear.verbose: false
 fem.solver.verbose: false
 fem.solver.newton.maxlineariterations: 100000
 fem.solver.newton.tolerance: 1e-07
-fem.solver.newton.linabstol: 1.25e-9
-fem.solver.newton.linreduction: 1.25e-9
+fem.solver.newton.linabstol: 2e-8
+fem.solver.newton.linreduction: 2e-8
+
 
 
 dgdiffusionflux.method: CDG2 # diffusion flux: CDG2, CDG, BR2, IP, NIPG, BO
-dgdiffusionflux.theoryparameters: 0 # scaling with theory parameters
-dgdiffusionflux.penalty: 20.
+dgdiffusionflux.theoryparameters: 1 # scaling with theory parameters
+dgdiffusionflux.penalty: 0.
 dgdiffusionflux.liftfactor: 1.0

pydemo/scalar/scalar.py

@@ -6,10 +6,7 @@ class Transport1D:
     def F_c(t,x,U):
         return as_matrix( [ [U[0], 0] ] )
 
-    def outflowValue(t,x,u):
-        return u
-    boundaryValue = {}
-    for i in range(1,5): boundaryValue.update( {i: outflowValue} )
+    boundary = {range(1,5): lambda t,x,u: u}
 
     def maxLambda(t,x,U,n):
         return abs(n[0])
@@ -48,19 +45,21 @@ def LinearAdvectionDiffusion1DMixed(v,eps,bnd):
     class Model(LinearAdvectionDiffusion1D(v,eps)):
         def dirichletValue(t,x,u):
             return bnd(t,x)
-        boundaryValue = {1:dirichletValue, 2:dirichletValue}
         def zeroFlux(t,x,u,n):
             return 0
-        diffusionBoundaryFlux = {3: zeroFlux, 4: zeroFlux}
+        if v is not None and eps is not None:
+            boundary = {(1,2): dirichletValue, (3,4): [zeroFlux,zeroFlux] }
+        else:
+            boundary = {(1,2): dirichletValue, (3,4): zeroFlux }
     return Model
 def LinearAdvectionDiffusion1DDirichlet(v,eps,bnd):
     class Model(LinearAdvectionDiffusion1D(v,eps)):
         def dirichletValue(t,x,u):
             return bnd(t,x)
-        boundaryValue = {}
-        for i in range(1,5): boundaryValue.update( {i: dirichletValue} )
+        boundary = { range(1,5): dirichletValue }
     return Model
 
+# burgers problems still need to be adapted to new API
 class Burgers1D:
     name = 'burgers'
 
@@ -96,13 +95,13 @@ def shockTransport():
            "shockTransport", None
 
 def sinProblem():
-    u0 = lambda t,x: as_vector( [sin(2*pi*x[0])] )
-    return LinearAdvectionDiffusion1DMixed(None,0.5,u0), u0(0,x),\
+    eps = 0.5
+    u0 = lambda t,x: as_vector( [sin(2*pi*x[0])*exp(-t*eps*(2*pi)**2)] )
+    return LinearAdvectionDiffusion1DMixed(None,eps,u0), u0(0,x),\
            [-1, 0], [1, 0.1], [50, 7], 0.2,\
-           "cos", None
+           "cos", lambda t: u0(t,x)
 
-def pulse():
-    eps = 0.001
+def pulse(eps=None):
     center = as_vector([ 0.5,0.5 ])
     x0 = x[0] - center[0]
     x1 = x[1] - center[1]
@@ -112,13 +111,18 @@ def pulse():
 
     def u0(t,x):
         sig2 = 0.004
-        sig2PlusDt4 = sig2+(4.0*eps*t)
+        if eps is None:
+            sig2PlusDt4 = sig2
+        else:
+            sig2PlusDt4 = sig2+(4.0*eps*t)
         xq = ( x0*cos(4.0*t) + x1*sin(4.0*t)) - 0.25
         yq = (-x0*sin(4.0*t) + x1*cos(4.0*t))
         return as_vector( [(sig2/ (sig2PlusDt4) ) * exp (-( xq*xq + yq*yq ) / sig2PlusDt4 )] )
     return LinearAdvectionDiffusion1DDirichlet([ux,uy],eps,u0), u0(0,x),\
            [0, 0], [1, 1], [16,16], 1.0,\
-           "pulse", lambda t: u0(t,x)
+           "pulse"+str(eps), lambda t: u0(t,x)
+def diffusivePulse():
+    return pulse(0.001)
 
 def burgersShock():
     return Burgers1D, riemanProblem(x[0],-0.5,[1],[0])

pydemo/scalar/shock_tube.py

@@ -8,9 +8,9 @@ from dune.femdg import createFemDGSolver
 from ufl import dot
 
 # from scalar import shockTransport as problem
-# from scalar import constantTransport as probelm
 # from scalar import sinProblem as problem
-from scalar import pulse as problem
+# from scalar import pulse as problem
+from scalar import diffusivePulse as problem
 
 Model, initial, x0,x1,N, endTime, name, exact = problem()
 

python/dune/femdg/_operators.py

@@ -273,13 +273,42 @@ def createFemDGSolver(Model, space,
             targs=['class Intersection, class Point'], static=True,
             predefined=predefined)
 
-    boundaryFluxDict = getattr(Model,"boundaryFlux",None)
-    if boundaryFluxDict is not None:
-        boundaryFlux = {}
-        for id,f in boundaryFluxDict.items(): boundaryFlux.update({id:f(t,x,u,n)})
-    else:
-        boundaryFlux = {}
-    generateMethod(struct, boundaryFlux,
+
+    #####################
+    ## boundary treatment
+    boundaryDict = getattr(Model,"boundary",{})
+    boundaryAFlux = {}
+    boundaryDFlux = {}
+    boundaryValue = {}
+    for k,f in boundaryDict.items():
+        # collect all ids (could be list or range)
+        ids = []
+        try:
+            for kk in k:
+                ids += [kk]
+        except TypeError:
+            ids += [k]
+        # TODO: check that id is not used more then once
+        # figure out what type of boundary condition is used
+        if isinstance(f,tuple) or isinstance(f,list):
+            assert hasAdvFlux and hasDiffFlux, "two boundary fluxes provided for id "+str(k)+" but only one bulk flux given"
+            method = [f[0](t,x,u,n), f[1](t,x,u,n)]
+            boundaryAFlux.update( [ (kk,method[0]) for kk in ids] )
+            boundaryDFlux.update( [ (kk,method[1]) for kk in ids] )
+        else:
+            try:
+                method = f(t,x,u,n)
+                if hasAdvFlux and not hasDiffFlux:
+                    boundaryAFlux.update( [ (kk,method) for kk in ids] )
+                elif not hasAdvFlux and hasDiffFlux:
+                    boundaryDFlux.update( [ (kk,method) for kk in ids] )
+                else:
+                    assert not (hasAdvFlux and hasDiffFlux), "one boundary fluxes provided for id "+str(k)+" but two bulk flux given"
+            except TypeError:
+                method = f(t,x,u)
+                boundaryValue.update( [ (kk,method) for kk in ids] )
+
+    generateMethod(struct, boundaryAFlux,
             'bool', 'boundaryFlux',
             args=['const int bndId',
                   'const double &t',
@@ -289,13 +318,7 @@ def createFemDGSolver(Model, space,
                   'RangeType &result'],
             targs=['class Entity, class Point'], static=True,
             predefined=predefined)
-    diffusionBoundaryFluxDict = getattr(Model,"diffusionBoundaryFlux",None)
-    if diffusionBoundaryFluxDict is not None:
-        diffusionBoundaryFlux = {}
-        for id,f in diffusionBoundaryFluxDict.items(): diffusionBoundaryFlux.update({id:f(t,x,u,n)})
-    else:
-        diffusionBoundaryFlux = {}
-    generateMethod(struct, diffusionBoundaryFlux,
+    generateMethod(struct, boundaryDFlux,
             'bool', 'diffusionBoundaryFlux',
             args=['const int bndId',
                   'const double &t',
@@ -306,12 +329,6 @@ def createFemDGSolver(Model, space,
                   'RangeType &result'],
             targs=['class Entity, class Point'], static=True,
             predefined=predefined)
-    boundaryValueDict = getattr(Model,"boundaryValue",None)
-    if boundaryValueDict is not None:
-        boundaryValue = {}
-        for id,f in boundaryValueDict.items(): boundaryValue.update({id:f(t,x,u)})
-    else:
-        boundaryValue = {}
     generateMethod(struct, boundaryValue,
             'bool', 'boundaryValue',
             args=['const int bndId',