extended implementation of dg stepper export - not quite working yet

dune/fem-dg/models/modelwrapper.hh

@@ -145,10 +145,12 @@ namespace Fem
     using BaseType::time_;
 
    public:
+#if 0
     ModelWrapper()
       : impl()
     {
     }
+#endif
 
     ModelWrapper( const ProblemType& problem )
       : impl_( problem )

dune/fem-dg/solver/dg.hh

@@ -26,7 +26,7 @@ namespace Fem
   template < class DestinationImp,
              class AdvectionModel,
              class DiffusionModel >
-  class DGOperator : public Fem::SpaceOperatorInterface< typename Traits::DestinationType >
+  class DGOperator : public Fem::SpaceOperatorInterface< DestinationImp >
   {
   public:
     typedef DestinationImp   DestinationType;
@@ -35,7 +35,7 @@ namespace Fem
     static const int polynomialOrder = DiscreteFunctionSpaceType :: polynomialOrder ;
 
     typedef typename DiscreteFunctionSpaceType :: GridPartType   GridPartType;
-    typedef typename GridPartType                                GridType;
+    typedef typename GridPartType::GridType                      GridType;
 
     typedef ModelImplementationWrapper< AdvectionModel >         ProblemType;
     typedef ModelWrapper< GridType, ProblemType >                ModelType;
@@ -55,10 +55,13 @@ namespace Fem
 
     static std::string name() const { return std::string("DGOperator"); }
 
-    DGOperator( const DiscreteFunctionSpaceType& space )
+    DGOperator( const DiscreteFunctionSpaceType& space,
+                const AdvectionModel &advectionModel,
+                const DiffusionModel &diffusionModel
+        )
       : extra_(),
         tp_(),
-        model_(),
+        model_(advectionModel),
         dgOperator( space.gridPart(), model_, extra_, name() ),
         rkSolver_( tp_, dgOperator_, dgOperator_, dgOperator_, name() )
     {}

pydemo/euler/euler.ufl

@@ -18,29 +18,16 @@ class Compressible2DEuler:
     def alpha(U,n):
         rho, u,v, p = Compressible2DEuler.toPrim(U)
         return abs(u*n[0]+v*n[1]) + sqrt(gamma*p/rho)
-class LLFFlux:
-    def flux(Descr,U,n,dt):
-        return jump(U)
-        UL = U('-')
-        UR = U('+')
-        return UL-UR
-        FL = Descr.F_c(UL)*n
-        FR = Descr.F_c(UR)*n
-        flux = (FL+FR)
-        max_value = lambda a,b: (a+b+abs(a-b))/2.
-        visc = max_value( Descr.alpha(UL,n), Descr.alpha(UR,n) )
-        return flux + visc/dt*(UR-UL)
 
 Model = Compressible2DEuler
-NumFLux = LLFFlux
 
+from dune.ufl import Space, NamedConstant
 space = Space(2,4)
 u = TrialFunction(space)
 v = TestFunction(space)
-dt = NamedConstant(triangle, "dt")    # time step
 
 # just a starting point...
 def femdgModel(description):
     return inner(description.F_c(u),grad(v))*dx
 
-F = -dt*femdgModel(Model)
+F = femdgModel(Model)

pydemo/euler/llf.py

+class LLFFlux:
+    def flux(Descr,U,n,dt):
+        UL = U('-')
+        UR = U('+')
+        FL = Descr.F_c(UL)*n('+')
+        FR = Descr.F_c(UR)*n('+')
+        flux = (FL+FR)/2
+        max_value = lambda a,b: (a+b+abs(a-b))/2.
+        visc = max_value( Descr.alpha(UL,n('-')), Descr.alpha(UR,n('+')) )
+        return flux + visc*(UR-UL)/2
+NumFlux = LLFFlux.flux

pydemo/euler/shock_tube.py

@@ -2,62 +2,80 @@ import time
 from dune.grid import structuredGrid
 from dune.fem import parameter
 import dune.create as create
-from ufl import TestFunction, TrialFunction, SpatialCoordinate, triangle, exp,\
-                FacetNormal, dx, grad, inner, as_vector, replace, sqrt,\
-                conditional, as_matrix, Max, dS, ds, avg, jump, outer
-from dune.ufl import NamedConstant
+from dune.models.elliptic.formfiles import loadModels
+from llf import NumFlux
+from ufl import *
 from dune.femdg import createFemDGSolver
 
+exec(open("euler.ufl").read())
+
 parameter.append({"fem.verboserank": -1})
 
 grid = structuredGrid([-1, 0], [1, 0.1], [200, 10])
-
-order     = 0
 dimR      = 4
-spaceName = "dgonb"
-space = create.space(spaceName, grid, dimrange=dimR, order=order)
+UL = as_vector( Model.toCons([1,0,0,1]) )
+UR = as_vector( Model.toCons([0.125,0,0,0.1]) )
+x = SpatialCoordinate(space.cell())
+initial = conditional(x[0]<0,UL,UR)
+t = 0
+saveStep = 0.01
+saveTime = saveStep
+count = 0
+endTime = 0.4
 
-import euler
+def useGalerkinOp():
+    # a full fv implementation using UFL and the galerkin operator
+    # some sign error or something else leading to wrong results
+    # still needs fixing
+    spaceName = "finitevolume"
+    space = create.space(spaceName, grid, dimrange=dimR)
 
-Euler = euler.Model
-NumFLux = euler.NumFlux
+    n = FacetNormal(space.cell())
 
-u = TrialFunction(space)
-v = TestFunction(space)
-dt = NamedConstant(triangle, "dt")    # time step
-x = SpatialCoordinate(space.cell())
-n = FacetNormal(space.cell())
+    u_h   = space.interpolate(initial, name='u_h')
+    u_h_n = u_h.copy(name="previous")
 
-UL = as_vector( Euler.toCons([1,0,0,1]) )
-UR = as_vector( Euler.toCons([0.125,0,0,0.1]) )
-initial = conditional(x[0]<0,UL,UR)
+    fullModel = inner( Model.F_c(u), grad(v) ) * dx -\
+                inner( NumFlux(Model,u,n,dt), v('+')-v('-')) * dS -\
+                inner( Model.F_c(u)*n, v) * ds
+    operator = create.operator("galerkin",
+                 inner(u-u_h_n,v)*dx == dt*replace(fullModel,{u:u_h_n}),
+                 space )
 
-u_h = space.interpolate(initial, name='u_h')
-u_h_n = u_h.copy(name="previous")
+    operator.model.dt = 1e-5
 
-# here a full implementation using UFL and the galerkin operator
-fullModel = inner(u_h_n,v)*dx-euler.F +\
-            dt*inner(LLFFlux.flux(Euler,u,n,dt),avg(v))*dS +\
-            dt*inner(Euler.F_c(u)*n,v)*ds
-# operator = create.operator("galerkin", fullModel, space )
-operator = createFemDGSolver( "euler", space, fullModel )
+    start = time.time()
+    grid.writeVTK('sod', pointdata=[u_h], number=count)
+    while t < 0.4:
+        u_h_n.assign(u_h)
+        operator.solve(target=u_h)
+        t += operator.model.dt
+        if t > saveTime:
+            count += 1
+            grid.writeVTK('sod', pointdata=[u_h], number=count)
+            saveTime += saveStep
+    grid.writeVTK('sod', pointdata=[u_h], number=count)
+    print("time loop:",time.time()-start)
 
-operator.model.dt = 5e-4
+def useODESolver():
+    spaceName = "dgonb"
+    polOrder = 0
+    space = create.space(spaceName, grid, order=polOrder, dimrange=dimR)
+    u_h   = space.interpolate(initial, name='u_h')
+    u_h_n = u_h.copy(name="previous")
+    operator = createFemDGSolver( Compressible2DEuler, space )
 
-t = 0
-saveStep = 0.01
-saveTime = saveStep
-count = 0
-grid.writeVTK('sod', pointdata=[u_h], number=count)
-
-start = time.time()
-while t < 0.4:
-    u_h_n.assign(u_h)
-    operator(u_h_n,u_h)
-    t += operator.model.dt
-    if t > saveTime:
-        count += 1
-        grid.writeVTK('sod', pointdata=[u_h], number=count)
-        saveTime += saveStep
-grid.writeVTK('sod', pointdata=[u_h], number=count)
-print("time loop:",time.time()-start)
+    start = time.time()
+    grid.writeVTK('sod', pointdata=[u_h], number=count)
+    while t < 0.4:
+        u_h_n.assign(u_h)
+        operator(u_h_n, u_h)
+        t += dt
+        if t > saveTime:
+            count += 1
+            grid.writeVTK('sod', pointdata=[u_h], number=count)
+            saveTime += saveStep
+    grid.writeVTK('sod', pointdata=[u_h], number=count)
+    print("time loop:",time.time()-start)
+
+useODESolver()

python/dune/femdg/_operators.py

@@ -21,8 +21,8 @@ def createLimiter(domainSpace, rangeSpace=None, bounds = [1e-12,1.], limiter='sc
     domainSpaceType = domainSpace._typeName
     rangeSpaceType = rangeSpace._typeName
 
-    _, domainFunctionIncludes, domainFunctionType, _, _ = domainSpace.storage
-    _, rangeFunctionIncludes, rangeFunctionType, _, _ = rangeSpace.storage
+    _, domainFunctionIncludes, domainFunctionType, _, _, _ = domainSpace.storage
+    _, rangeFunctionIncludes, rangeFunctionType, _, _, _ = rangeSpace.storage
 
     includes = ["dune/fem-dg/operator/limiter/limiter.hh"]
     includes += domainSpace._includes + domainFunctionIncludes
@@ -70,30 +70,57 @@ def createOrderRedcution(domainSpace):
 
 
 # create DG operator + solver
-def createFemDGSolver(name, space, advectionModel, diffusionModel = None ):
-
-    if diffusionModel is None:
-        diffusionModel = advectionModel
+def createFemDGSolver(Model, space):
+    from ufl import TestFunction,TrialFunction,dx,grad,inner,zero
+    import dune.create as create
+
+    u = TrialFunction(space)
+    v = TestFunction(space)
+    if hasattr(Model,"F_c"):
+        advModel = inner(Model.F_c(u),grad(v))*dx
+    else:
+        advModel = inner(grad(u-u),grad(v))*dx
+    if hasattr(Model,"S_ns"):
+        advModel += inner(S_ns(u,grad(u)),v)*dx
+    if hasattr(Model,"F_d"):
+        diffModel = inner(Model.F_d(u,grad(u)),grad(v))*dx
+    else:
+        diffModel = inner(grad(u-u),grad(v))*dx
+    if hasattr(Model,"S_s"):
+        diffModel += inner(S_s(u,grad(u)),v)*dx
+    advModel  = create.model("elliptic",space.grid, advModel)
+    diffModel = create.model("elliptic",space.grid, diffModel)
 
     spaceType = space._typeName
 
-    advectionModelType = advectionModel._typeName
-    diffusionModelType = diffusionModel._typeName
+    modelType = "DiffusionModel< " +\
+          "typename " + spaceType + "::GridPartType, " +\
+          spaceType + "::dimRange, " +\
+          spaceType + "::dimRange, " +\
+          "typename " + spaceType + "::RangeFieldType >"
 
-    _, destinationIncludes, destinationType, _, _ = space.storage
+    advModelType  = modelType
+    diffModelType = modelType
 
-    print("model name :", advectionModelType )
+    _, destinationIncludes, destinationType, _, _, _ = space.storage
 
-    includes  = [ name + '.hh' ]
-    includes += ["dune/fem-dg/solver/dg.hh"]
+    includes  = ["dune/fem-dg/solver/dg.hh"]
     includes += space._includes + destinationIncludes
+    includes += ["dune/fem/schemes/diffusionmodel.hh", "dune/fempy/parameter.hh"]
 
-    typeName = 'Dune::Fem::DGOperator< ' + destinationType + ', ' + advectionModelType + ', ' + diffusionModelType + ' >'
+    typeName = 'Dune::Fem::DGOperator< ' + destinationType + ', ' + advModelType + ', ' + diffModelType + ' >'
 
-    constructor = Constructor(['const '+spaceType + ' &space'],
+    constructor = Constructor(['const '+spaceType + ' &space',
+                               'const '+advModelType + ' &advectionModel',
+                               'const '+diffModelType + ' &diffusionModel'
+                              ],
                               ['return new ' + typeName + '(space);'],
                               ['"space"_a',
-                               'pybind11::keep_alive< 1, 2 >()', 'pybind11::keep_alive< 1, 3 >()'])
+                               '"advectionModel"_a',
+                               '"diffusionModel"_a',
+                               'pybind11::keep_alive< 1, 2 >()',
+                               'pybind11::keep_alive< 1, 3 >()',
+                               'pybind11::keep_alive< 1, 4 >()'])
 
     # add method activated to inspect limited cells.
     setTimeStepSize = Method('setTimeStepSize', '&'+typeName+'::setTimeStepSize')