make scalar advection diffusion problems work

mainly: added diffusion boundary flux. The boundary treatment needs some
overhaul

dune/fem-dg/models/modelwrapper.hh

@@ -217,6 +217,7 @@ namespace Fem
                            JacobianRangeType& f ) const
     {
       assert( hasAdvection );
+      advection_.init( local.entity() ); // TODO: this should not be required but the pulse problem fails without
       advection_.diffusiveFlux( local.quadraturePoint(), u, du, f);
     }
 
@@ -246,6 +247,7 @@ namespace Fem
                            RangeType& A ) const
     {
       assert( hasAdvection );
+      assert( 0 ); // if this is used we have to check if this is correct
       // TODO: u != ubar and du != dubar
       advection_.linDiffusiveFlux( u, du, local.quadraturePoint(), u, du, A);
     }
@@ -292,12 +294,6 @@ namespace Fem
       const bool isFluxBnd =
 #endif
       AdditionalType::boundaryFlux(id, time(), local.entity(), local.quadraturePoint(), normal, uLeft, gLeft);
-#if 0
-      std::cout << "BoundaryFlux with id=" << id << " with uLeft=" << uLeft
-                << " results in flux=" << gLeft
-                << " --- " << isFluxBnd
-                << std::endl;
-#endif
       assert( isFluxBnd );
       return 0; // QUESTION: do something better here? Yes, return time step restriction if possible
     }
@@ -321,9 +317,14 @@ namespace Fem
                                          const JacobianRangeType& jacLeft,
                                          RangeType& gLeft ) const
     {
-      // TODO: need to add something to 'Addtional'?
-      assert( hasDiffusion );
-      return 0;
+      const DomainType normal = local.intersection().integrationOuterNormal( local.localPosition() );
+      int id = getBoundaryId( local );
+#ifndef NDEBUG
+      const bool isFluxBnd =
+#endif
+      AdditionalType::diffusionBoundaryFlux(id, time(), local.entity(), local.quadraturePoint(), normal, uLeft, jacLeft, gLeft);
+      assert( isFluxBnd );
+      return 0; // QUESTION: do something better here? Yes, return time step restriction if possible
     }
 
     template <class LocalEvaluation>

dune/fem-dg/operator/dg/passtraits.hh

@@ -57,13 +57,16 @@ namespace Fem
     typedef Fem::CachingQuadrature< GridPartType, 1 >    FaceQuadratureType;
 
     // Allow generalization to systems
+#if 1
     typedef Fem::DiscontinuousGalerkinSpace<
                                         FunctionSpaceType,
                                         GridPartType, polynomialOrder,
                                         Fem::CachingStorage >             DiscreteFunctionSpaceType;
-    //typedef Fem::LegendreDiscontinuousGalerkinSpace< FunctionSpaceType,
-    //                                    GridPartType, polOrd,
-    //                                    Fem::CachingStorage >             DiscreteFunctionSpaceType;
+#else
+    typedef Fem::LegendreDiscontinuousGalerkinSpace< FunctionSpaceType,
+                                        GridPartType, polOrd,
+                                        Fem::CachingStorage >             DiscreteFunctionSpaceType;
+#endif
     typedef Fem::AdaptiveDiscreteFunction< DiscreteFunctionSpaceType >    DestinationType;
 
     // Indicator for Limiter

pydemo/euler/shock_tube.py

@@ -73,9 +73,12 @@ def useGalerkinOp():
 
 def useODESolver(polOrder=2, limiter='default'):
     global count, t, dt, saveTime
-    spaceName = "dgonb"
     polOrder = polOrder
-    space = create.space(spaceName, grid, order=polOrder, dimrange=dimR)
+    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 )
@@ -98,7 +101,7 @@ def useODESolver(polOrder=2, limiter='default'):
         t += dt
         tcount += 1
         if t > saveTime:
-            print('dt = ', dt, 'time = ',t, 'count = ',count )
+            print('dt = ', dt, 'time = ',t, 'count = ',count, flush=True )
             count += 1
             # rho, v, p = Model.toPrim(u_h) # is this needed - works for me # without and it should...
             grid.writeVTK(name,
@@ -120,7 +123,6 @@ if True:
     grid = structuredGrid(x0,x1,N)
     # grid = create.grid("ALUSimplex", cartesianDomain(x0,x1,N))
     useODESolver(2,'default')      # third order with limiter
-#    useODESolver(2,None)      # third order with limiter
 elif False:
     N = [n*10 for n in N]
     grid = structuredGrid(x0,x1,N)

pydemo/scalar/paramSolver

@@ -13,18 +13,18 @@ fem.timeprovider.factor: 0.45
 fem.timeprovider.updatestep: 1
 
 # parameter for the implicit solvers
-fem.differenceoperator.eps: 1e-12
+# fem.differenceoperator.eps: 1e-12
 fem.solver.gmres.restart: 50
 fem.solver.newton.verbose: true
-fem.solver.newton.linear.verbose: false
-fem.solver.verbose: false
-fem.solver.newton.maxlineariterations: 1000
-fem.solver.newton.tolerance: 1e-10
-fem.solver.newton.linabstol: 1.25e-12
-fem.solver.newton.linreduction: 1.25e-12
+fem.solver.newton.linear.verbose: true
+fem.solver.verbose: true
+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
 
 
 dgdiffusionflux.method: CDG2 # diffusion flux: CDG2, CDG, BR2, IP, NIPG, BO
-dgdiffusionflux.theoryparameters: 1 # scaling with theory parameters
-dgdiffusionflux.penalty: 0.
+dgdiffusionflux.theoryparameters: 0 # scaling with theory parameters
+dgdiffusionflux.penalty: 20.
 dgdiffusionflux.liftfactor: 1.0

pydemo/scalar/scalar.py

@@ -20,33 +20,46 @@ class Transport1D:
     def jump(U,V):
         return U-V
 
-class LinearAdvectionDiffusion1D:
-    name = 'linearAD'
-
-    # def F_c(t,x,U):
-    #     return as_matrix( [ [U[0], 0] ] )
-    def F_v(t,x,U,DU):
-        return 0.5*DU
-    def maxLambda(t,x,U,n):
-        return abs(n[0])
-    def velocity(t,x,U):
-        return as_vector( [0,0] )
-    def physical(U):
-        return 1
-    def jump(U,V):
-        return U-V
-def LinearAdvectionDiffusion1DDirichlet(bnd):
-    class Model(LinearAdvectionDiffusion1D):
+def LinearAdvectionDiffusion1D(v,eps):
+    v = as_vector(v)
+    class Model:
+        if v is not None and eps is not None:
+            name = 'linearAD'
+        elif v is not None:
+            name = 'linearA'
+        else:
+            name = 'linearD'
+        if v is not None:
+            def F_c(t,x,U):
+                return as_matrix( [[ *(v*U[0]) ]] )
+        if eps is not None:
+            def F_v(t,x,U,DU):
+                return eps*DU
+        def maxLambda(t,x,U,n):
+            return abs(dot(v,n))
+        def velocity(t,x,U):
+            return v
+        def physical(U):
+            return 1
+        def jump(U,V):
+            return U-V
+    return Model
+def LinearAdvectionDiffusion1DMixed(v,eps,bnd):
+    class Model(LinearAdvectionDiffusion1D(v,eps)):
+        def dirichletValue(t,x,u):
+            return bnd
+        boundaryValue = {1:dirichletValue, 2:dirichletValue}
+        def zeroFlux(t,x,u,n):
+            return 0
+        diffusionBoundaryFlux = {3: zeroFlux, 4: zeroFlux}
+    return Model
+def LinearAdvectionDiffusion1DDirichlet(v,eps,bnd):
+    class Model(LinearAdvectionDiffusion1D(v,eps)):
         def dirichletValue(t,x,u):
             return bnd
         boundaryValue = {}
         for i in range(1,5): boundaryValue.update( {i: dirichletValue} )
     return Model
-class LinearAdvectionDiffusion1DNeuman(LinearAdvectionDiffusion1D):
-    def zeroFlux(t,x,u,n):
-        return 0
-    boundaryFlux = {}
-    for i in range(1,5): boundaryFlux.update( {i: zeroFlux} )
 
 class Burgers1D:
     name = 'burgers'
@@ -74,13 +87,19 @@ def riemanProblem(x,x0,UL,UR):
     return conditional(x<x0,as_vector(UL),as_vector(UR))
 
 def constantTransport():
-    return Transport1D, as_vector( [0.1] )
+    return Transport1D, as_vector( [0.1] ),\
+           [-1, 0], [1, 0.1], [50, 7], 0.1,\
+           "constant"
 def shockTransport():
-    return Transport1D, riemanProblem(x[0],-0.5, [1], [0])
+    return Transport1D, riemanProblem(x[0],-0.5, [1], [0]),\
+           [-1, 0], [1, 0.1], [50, 7], 1.0,\
+           "shockTransport"
 
 def sinProblem():
     u0 = as_vector( [sin(2*pi*x[0])] )
-    return LinearAdvectionDiffusion1DDirichlet(u0), u0
+    return LinearAdvectionDiffusion1DMixed(None,0.5,u0), u0,\
+           [-1, 0], [1, 0.1], [50, 7], 0.2,\
+           "cos"
 
 def pulse():
     t = 0
@@ -89,15 +108,17 @@ def pulse():
     x0 = x[0] - center[0]
     x1 = x[1] - center[1]
     sig2 = 0.004
-    sig2PlusDt4 = sig2+(4.0*eps*t);
-    xq = ( x0*cos(4.0*t) + x1*sin(4.0*t))
+    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))
 
-    u0 = as_vector( [(sig2/ (sig2PlusDt4) ) * exp (-( xq*xq + yq*yq ) / sig2PlusDt4 )] )
-    return LinearAdvectionDiffusion1DDirichlet(u0), u0
+    ux = -4.0*x1;
+    uy =  4.0*x0;
 
-def cosProblem():
-    return LinearAdvectionDiffusion1DNeuman, as_vector( [cos(2*pi*x[0])] )
+    u0 = as_vector( [(sig2/ (sig2PlusDt4) ) * exp (-( xq*xq + yq*yq ) / sig2PlusDt4 )] )
+    return LinearAdvectionDiffusion1DDirichlet([ux,uy],eps,u0), u0,\
+           [0, 0], [1, 1], [8,8], 1.0,\
+           "puls"
 
 def burgersShock():
     return Burgers1D, riemanProblem(x[0],-0.5,[1],[0])

pydemo/scalar/shock_tube.py

@@ -7,19 +7,14 @@ from dune.femdg import createFemDGSolver
 
 # from scalar import shockTransport as problem
 # from scalar import constantTransport as probelm
-from scalar import sinProblem as problem
-# from scalar import cosProblem as problem
-# from scalar import burgersShock as problem
-# from scalar import burgersVW as problem
-# from scalar import burgersStationaryShock as problem
+# from scalar import sinProblem as problem
+from scalar import pulse as problem
 
-Model, initial = problem()
+Model, initial, x0,x1,N, endTime, name = problem()
 
 parameter.append({"fem.verboserank": 0})
 parameter.append("parameter")
 
-x0,x1,N = [-1, 0], [1, 0.1], [50, 7]
-#x0,x1,N = [0, 0], [1, 1], [8,8]
 grid = structuredGrid(x0,x1,N)
 
 # grid = create.grid("ALUSimplex", cartesianDomain(x0,x1,N))
@@ -29,7 +24,6 @@ dt       = 0.001
 saveStep = 0.01
 saveTime = saveStep
 count    = 0
-endTime  = 1.
 
 def useODESolver():
     global count, t, dt, saveTime
@@ -38,22 +32,22 @@ def useODESolver():
     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 )
+    operator = createFemDGSolver( Model, space, limiter=None, diffusionScheme="cdg2" )
 
     operator.setTimeStepSize(dt)
 
     start = time.time()
-    grid.writeVTK(Model.name, celldata=[u_h], number=count)
+    grid.writeVTK(name, celldata=[u_h], number=count, subsampling=2)
     while t < endTime:
         operator.solve(target=u_h)
         dt = operator.deltaT()
         t += dt
         if t > saveTime:
-            print('dt = ', dt, 'time = ',t, 'count = ',count )
+            print('dt = ', dt, 'time = ',t, 'count = ',count, flush=True )
             count += 1
-            grid.writeVTK(Model.name, celldata=[u_h], number=count)
+            grid.writeVTK(name, celldata=[u_h], number=count, subsampling=2)
             saveTime += saveStep
-    grid.writeVTK(Model.name, celldata=[u_h], number=count)
+    grid.writeVTK(name, celldata=[u_h], number=count, subsampling=2)
     print("time loop:",time.time()-start)
 
 useODESolver()

python/dune/femdg/_operators.py

@@ -149,7 +149,7 @@ def generateMethod(struct,expr, cppType, name,
 
 # create DG operator + solver
 def createFemDGSolver(Model, space,
-        limiter="default"):
+        limiter="default", diffusionScheme = "cdg2"):
     import dune.create as create
 
     u = TrialFunction(space)
@@ -166,8 +166,6 @@ def createFemDGSolver(Model, space,
     hasNonStiffSource = hasattr(Model,"S_ns")
     if hasNonStiffSource:
         advModel += inner(as_vector(S_ns(t,x,u,grad(u))),v)*dx
-    else:
-        advModel += inner(t*u,v)*dx
 
     hasDiffFlux = hasattr(Model,"F_v")
     if hasDiffFlux:
@@ -177,8 +175,6 @@ def createFemDGSolver(Model, space,
     hasStiffSource = hasattr(Model,"S_s")
     if hasStiffSource:
         diffModel += inner(as_vector(S_s(t,x,u,grad(u))),v)*dx
-    else:
-        advModel += inner(t*u,v)*dx
 
     advModel  = create.model("elliptic",space.grid, advModel)
     diffModel = create.model("elliptic",space.grid, diffModel)
@@ -270,6 +266,23 @@ 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,
+            'bool', 'diffusionBoundaryFlux',
+            args=['const int bndId',
+                  'const double &t',
+                  'const Entity& entity', 'const Point &x',
+                  'const DomainType &normal',
+                  'const RangeType &u',
+                  'const JacobianRangeType &jac',
+                  'RangeType &result'],
+            targs=['class Entity, class Point'], static=True,
+            predefined=predefined)
     boundaryValueDict = getattr(Model,"boundaryValue",None)
     if boundaryValueDict is not None:
         boundaryValue = {}
@@ -330,12 +343,13 @@ def createFemDGSolver(Model, space,
     solverId   = "Dune::Fem::Solver::Enum::fem"
     formId     = "Dune::Fem::Formulation::Enum::primal"
     limiterId  = "Dune::Fem::AdvectionLimiter::Enum::limited"
-    advFluxId  = "Dune::Fem::AdvectionFlux::Enum::llf"
-
+    advFluxId  = "Dune::Fem::AdvectionFlux::Enum::none"
     diffFluxId = "Dune::Fem::DiffusionFlux::Enum::none"
-    if hasDiffFlux:
-        diffFluxId = "Dune::Fem::DiffusionFlux::Enum::primal"
 
+    if hasDiffFlux:
+        diffFluxId = "Dune::Fem::DiffusionFlux::Enum::"+diffusionScheme
+    if hasAdvFlux:
+        advFluxId  = "Dune::Fem::AdvectionFlux::Enum::llf"
     if limiter == None or limiter == False or limiter.lower() == "unlimiter":
         limiterId = "Dune::Fem::AdvectionLimiter::Enum::unlimited"