Make advection-diffusion work.

pydemo/scalar/paramSolver

@@ -5,17 +5,17 @@ fem.ode.odesolver: IM # 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: 5
-fem.ode.maxiterations: 100005
+fem.ode.miniterations: 35
+fem.ode.maxiterations: 100
 fem.ode.cflStart: 1.
-#fem.ode.cflMax: 5
+fem.ode.cflMax: 5
 fem.timeprovider.factor: 0.45
 fem.timeprovider.updatestep: 1
 
 # parameter for the implicit solvers
 # fem.differenceoperator.eps: 1e-12
 fem.solver.gmres.restart: 50
-fem.solver.newton.verbose: true
+fem.solver.newton.verbose: false
 fem.solver.newton.linear.verbose: false
 fem.solver.verbose: false
 fem.solver.newton.maxlineariterations: 100000

pydemo/scalar/scalar.py

@@ -36,10 +36,12 @@ 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
         def physical(U):
             return 1
         def jump(U,V):
-            return U-V
+            return abs(U-V)
     return Model
 def LinearAdvectionDiffusion1DMixed(v,eps,bnd):
     class Model(LinearAdvectionDiffusion1D(v,eps)):
@@ -102,7 +104,7 @@ def sinProblem():
            "cos", lambda t: u0(t,x)
 
 def pulse(eps=None):
-    center = as_vector([ 0.5,0.5 ])
+    center  = as_vector([ 0.5,0.5 ])
     x0 = x[0] - center[0]
     x1 = x[1] - center[1]
 
@@ -115,7 +117,7 @@ def pulse(eps=None):
             sig2PlusDt4 = sig2
         else:
             sig2PlusDt4 = sig2+(4.0*eps*t)
-        xq = ( x0*cos(4.0*t) + x1*sin(4.0*t)) - 0.25
+        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),\

pydemo/scalar/shock_tube.py

@@ -8,7 +8,7 @@ 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 pulse as problem
 from scalar import diffusivePulse as problem
 
@@ -22,8 +22,8 @@ grid = structuredGrid(x0,x1,N)
 # grid = create.grid("ALUSimplex", cartesianDomain(x0,x1,N))
 dimR     = 1
 t        = 0
-dt       = 0.005
-saveStep = 0.05
+dt       = 0.001
+saveStep = 0.1
 saveTime = saveStep
 count    = 0
 
@@ -36,20 +36,25 @@ def useODESolver():
     u_h   = space.interpolate(initial, name='u_h')
     operator = createFemDGSolver( Model, space, limiter=None, diffusionScheme="cdg2" )
 
-    operator.setTimeStepSize(dt)
+    #operator.setTimeStepSize(dt)
 
     start = time.time()
-    grid.writeVTK(name, celldata=[u_h], number=count, subsampling=2)
+    #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
         if t > saveTime:
-            print('dt = ', dt, 'time = ',t, 'count = ',count, flush=True )
             count += 1
-            grid.writeVTK(name, celldata=[u_h], number=count, subsampling=2)
+            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, subsampling=2)
+    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 )