fix an issue with the diffusive boundary flux (grad(u) wasn't used)

Highlight an issue with the diffusion timestep for IMEX method

pydemo/scalar/paramSolver

 # SOLVER CONFIGURATION
 #---------------------
 
-fem.ode.odesolver: IM # ode solvers: EX, IM, IMEX
+fem.ode.odesolver: IMEX # ode solvers: EX, IM, IMEX
 fem.ode.order: 3
 fem.ode.verbose: cfl # ode output: none, cfl, full
 fem.ode.cflincrease: 1.25

pydemo/scalar/scalar.py

@@ -36,8 +36,11 @@ def LinearAdvectionDiffusion1D(v,eps):
         if eps is not None:
             def F_v(t,x,U,DU):
                 return eps*DU
-            def maxDiffusion(t,x,U):
-                return eps
+            # TODO: this method is used in an IMEX method allough the
+            # diffusion is explicit - should we change that behaviour?
+            # commented out for test of IMEX
+            # def maxDiffusion(t,x,U):
+            #    return eps
         def physical(U):
             return 1
         def jump(U,V):
@@ -49,8 +52,10 @@ def LinearAdvectionDiffusion1DMixed(v,eps,bnd):
             return bnd(t,x)
         def zeroFlux(t,x,u,n):
             return 0
+        def zeroDiffFlux(t,x,u,du,n):
+            return 0
         if v is not None and eps is not None:
-            boundary = {(1,2): dirichletValue, (3,4): [zeroFlux,zeroFlux] }
+            boundary = {(1,2): dirichletValue, (3,4): [zeroFlux,zeroDiffFlux] }
         else:
             boundary = {(1,2): dirichletValue, (3,4): zeroFlux }
     return Model
@@ -101,7 +106,15 @@ def sinProblem():
     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", lambda t: u0(t,x)
+           "sin", lambda t: u0(t,x)
+
+def sinTransportProblem():
+    eps = 0.5
+    v   = [1,0]
+    u0 = lambda t,x: as_vector( [sin(2*pi*x[0])*exp(-t*eps*(2*pi)**2)] )
+    return LinearAdvectionDiffusion1DMixed(v,eps,u0), u0(0,x),\
+           [-1, 0], [1, 0.1], [50, 7], 0.2,\
+           "sin", lambda t: u0(t,x)
 
 def pulse(eps=None):
     center  = as_vector([ 0.5,0.5 ])

pydemo/scalar/shock_tube.py

@@ -8,9 +8,10 @@ from dune.femdg import createFemDGSolver
 from ufl import dot
 
 # from scalar import shockTransport as problem
-#from scalar import sinProblem 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
+# from scalar import diffusivePulse as problem
 
 Model, initial, x0,x1,N, endTime, name, exact = problem()
 
@@ -49,9 +50,9 @@ def useODESolver():
             exit(0)
         dt = operator.deltaT()
         t += dt
+        print('dt = ', dt, 'time = ',t, 'count = ',count, flush=True )
         if t > saveTime:
             count += 1
-            print('dt = ', dt, 'time = ',t, 'count = ',count, flush=True )
             grid.writeVTK(name, celldata=[u_h], number=count)
             saveTime += saveStep
     grid.writeVTK(name, celldata=[u_h], number=count)

python/dune/femdg/_operators.py

@@ -183,7 +183,7 @@ def createFemDGSolver(Model, space,
         advModel = inner(t*grad(u-u),grad(v))*dx    # TODO: make a better empty model
     hasNonStiffSource = hasattr(Model,"S_ns")
     if hasNonStiffSource:
-        advModel += inner(as_vector(S_ns(t,x,u,grad(u))),v)*dx
+        advModel += inner(as_vector(Model.S_ns(t,x,u,grad(u))),v)*dx
 
     hasDiffFlux = hasattr(Model,"F_v")
     if hasDiffFlux:
@@ -192,7 +192,7 @@ def createFemDGSolver(Model, space,
         diffModel = inner(t*grad(u-u),grad(v))*dx   # TODO: make a better empty model
     hasStiffSource = hasattr(Model,"S_s")
     if hasStiffSource:
-        diffModel += inner(as_vector(S_s(t,x,u,grad(u))),v)*dx
+        diffModel += inner(as_vector(Model.S_s(t,x,u,grad(u))),v)*dx
 
     advModel  = create.model("elliptic",space.grid, advModel)
     diffModel = create.model("elliptic",space.grid, diffModel)
@@ -329,7 +329,7 @@ def createFemDGSolver(Model, space,
         # 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)]
+            method = [f[0](t,x,u,n), f[1](t,x,u,grad(u),n)]
             boundaryAFlux.update( [ (kk,method[0]) for kk in ids] )
             boundaryDFlux.update( [ (kk,method[1]) for kk in ids] )
         else: