implemented some explicit RK schemes

dune/fem-dg/python/operator.hh

@@ -38,7 +38,9 @@ namespace Dune
            pybind11::keep_alive< 1, 5 >()
            );
       cls.def( "applyLimiter", []( Operator &self, DF &u) { self.applyLimiter(u); } );
-
+      cls.def( "setTime", &Operator::setTime);
+      cls.def_property_readonly("timeStepEstimate",
+          [](const Operator &self) -> double {return self.timeStepEstimate(); });
       Dune::Python::insertClass<ExplType>(cls,"ExplType",
           Dune::Python::GenerateTypeName("NotAvailable"),
           Dune::Python::IncludeFiles{});

dune/fem-dg/solver/rungekuttasolver.hh

@@ -342,11 +342,14 @@ namespace Fem
       // if it is not 0 otherwise use the internal estimate
       tpPtr_->provideTimeStepEstimate(maxTimeStep);
       // adjust fixed time step with timeprovider.factor()
+      // SAD: how is fixedTimeStep_ initialized here?
+      fixedTimeStep_ = 0;
       fixedTimeStep_ /= tpPtr_->factor() ;
       if ( fixedTimeStep_ > 1e-20 )
         tpPtr_->init( fixedTimeStep_ );
       else
         tpPtr_->init();
+      tpPtr_->init();
       std::cout << "cfl = " << double(tpPtr_->factor()) << ", T_0 = " << tpPtr_->time() << std::endl;
     }
 

pydemo/scalar/main.py

@@ -3,9 +3,9 @@ from dune.femdg.testing import run
 
 # from scalar import shockTransport as problem
 # from scalar import sinProblem as problem
-from scalar import sinTransportProblem as problem
+# from scalar import sinTransportProblem as problem
 # from scalar import sinAdvDiffProblem as problem
-# from scalar import pulse as problem
+from scalar import pulse as problem
 # from scalar import diffusivePulse as problem
 
 parameter.append({"fem.verboserank": 0})
@@ -26,21 +26,21 @@ parameters = {"fem.ode.odesolver": "EX",   # EX, IM, IMEX
               "dgdiffusionflux.liftfactor": 1}
 
 parameters['fem.ode.odesolver'] = 'EX'
-uh,errorEx = run(*problem(),
+uh,errorEx = run(problem(),
         startLevel=0, polOrder=2, limiter=None,
         primitive=None, saveStep=0.01, subsamp=0,
-        dt=0.001,
+        dt=None,
         parameters=parameters)
 parameters['fem.ode.odesolver'] = 'IM'
-uh,errorIm = run(*problem(),
+uh,errorIm = run(problem(),
         startLevel=0, polOrder=2, limiter=None,
         primitive=None, saveStep=0.01, subsamp=0,
-        dt=None, # 0.001, otherwise this is fixed time step
+        dt=None,
         parameters=parameters)
 parameters['fem.ode.odesolver'] = 'IMEX'
-uh,errorImex = run(*problem(),
+uh,errorImex = run(problem(),
         startLevel=0, polOrder=2, limiter=None,
         primitive=None, saveStep=0.01, subsamp=0,
-        dt=0.001,
+        dt=None,
         parameters=parameters)
 print(errorEx,errorIm,errorImex)

pydemo/scalar/rk.py

+import time, math, sys
+from dune.grid import structuredGrid
+from dune.fem.space import dgonb
+from dune.fem.function import integrate
+from dune.ufl import Constant
+from ufl import dot
+from dune.femdg import femDGOperator
+
+from dune.fem import parameter
+
+# from scalar import shockTransport as problem
+# from scalar import sinProblem as problem
+# from scalar import sinTransportProblem as problem
+# from scalar import sinAdvDiffProblem as problem
+from scalar import pulse as problem
+# from scalar import diffusivePulse as problem
+
+parameter.append({"fem.verboserank": 0})
+parameters = {"fem.ode.odesolver": "EX",   # EX, IM, IMEX
+              "fem.ode.order": 3,
+              "fem.ode.verbose": "cfl",      # none, cfl, full
+              "fem.ode.cflMax": 100,
+              "fem.ode.cflincrease": 1.25,
+              "fem.ode.miniterations": 35,
+              "fem.ode.maxiterations": 100,
+              "fem.timeprovider.factor": 0.45,
+              "dgadvectionflux.method": "LLF",
+              "fem.solver.gmres.restart": 50,
+              "dgdiffusionflux.method": "CDG2",      # CDG2, CDG, BR2, IP, NIPG, BO
+              "dgdiffusionflux.theoryparameters": 1, # scaling with theory parameters
+              "dgdiffusionflux.penalty": 0,
+              "dgdiffusionflux.liftfactor": 1}
+
+# http://www.sspsite.org
+# https://arxiv.org/pdf/1605.02429.pdf
+# https://openaccess.leidenuniv.nl/bitstream/handle/1887/3295/02.pdf?sequence=7
+# https://epubs.siam.org/doi/10.1137/07070485X
+
+class Heun:
+    def __init__(self, op, cfl=None):
+        self.stages = 2
+        self.op = op
+        self.c = [0,1]
+        self.b = [0.5,0.5]
+        self.A = [[0,0],[1,0]]
+        self.k = self.stages*[None]
+        self.cfl = 0.45 if cfl is None else cfl
+        for i in range(self.stages):
+            self.k[i] = op.space.interpolate(op.space.dimRange*[0],name="k")
+        self.tmp = op.space.interpolate(op.space.dimRange*[0],name="tmp")
+    def __call__(self,u,dt=None):
+        self.op(u,self.k[0])
+        if dt is None:
+            dt = self.op.timeStepEstimate*self.cfl
+        for i in range(1,self.stages):
+            self.tmp.assign(u)
+            for j in range(i):
+                self.tmp.axpy(dt*self.A[i][j],self.k[j])
+            self.op.stepTime(self.c[i],dt)
+            self.op(self.tmp,self.k[i])
+        for i in range(self.stages):
+            u.axpy(dt*self.b[i],self.k[i])
+        return dt
+
+class ExplSSP3:
+    def __init__(self, stages,op,cfl=None):
+        self.op     = op
+        self.n      = math.ceil(math.sqrt(stages))
+        self.stages = self.n*self.n
+        self.r      = self.stages-self.n
+        self.q2     = op.space.interpolate(space.dimRange*[0],name="q2")
+        self.tmp    = self.q2.copy()
+        self.cfl    = 0.45 * stages*(1-1/self.n) \
+                      if cfl is None else cfl
+        for i in range(1,self.stages+1):
+            print(i,self.c(i))
+    def c(self,i):
+        return (i-1)/(self.n*self.n-self.n) \
+               if i<=(self.n+2)*(self.n-1)/2 \
+               else (i-self.n)/(self.n*self.n-self.n)
+    def __call__(self,u,dt=None):
+        if dt is None:
+            self.op(u, self.tmp)
+            dt = self.op.timeStepEstimate*self.cfl
+        fac = dt/self.r
+        i = 1
+        while i <= (self.n-1)*(self.n-2)/2:
+            self.op.stepTime(self.c(i),dt)
+            self.op(u,self.tmp)
+            u.axpy(fac, self.tmp)
+            i += 1
+        self.q2.assign(u)
+        while i <= self.n*(self.n+1)/2:
+            self.op.stepTime(self.c(i),dt)
+            self.op(u,self.tmp)
+            u.axpy(fac, self.tmp)
+            i += 1
+        u.as_numpy[:] *= (self.n-1)/(2*self.n-1)
+        u.axpy(self.n/(2*self.n-1), self.q2)
+        while i <= self.stages:
+            self.op.stepTime(self.c(i),dt)
+            self.op(u,self.tmp)
+            u.axpy(fac, self.tmp)
+            i += 1
+        return dt
+class ExplSSP4_10:
+    def __init__(self, op,cfl=None):
+        self.op     = op
+        self.stages = 10
+        self.q2     = op.space.interpolate(space.dimRange*[0],name="q2")
+        self.tmp    = self.q2.copy()
+        self.cfl    = 0.45 * self.stages*0.6 if cfl is None else cfl
+        for i in range(1,self.stages+1):
+            print(i,self.c(i))
+    def c(self,i):
+        return (i-1)/6 if i<=5 else (i-4)/6
+    def __call__(self,u,dt=None):
+        if dt is None:
+            self.op(u, self.tmp)
+            dt = self.op.timeStepEstimate*self.cfl # how is this reset for the next time step
+        i = 1
+        self.q2.assign(u)
+        while i <= 5:
+            self.op.stepTime(self.c(i),dt)
+            self.op(u,self.tmp)
+            u.axpy(dt/6, self.tmp)
+            i += 1
+
+        self.q2.as_numpy[:] *= 1/25
+        self.q2.axpy(9/25,u)
+        u.as_numpy[:] *= -5
+        u.axpy(15,self.q2)
+
+        while i <= 9:
+            self.op.stepTime(self.c(i),dt)
+            self.op(u,self.tmp)
+            u.axpy(dt/6, self.tmp)
+            i += 1
+
+        self.op.stepTime(self.c(i),dt)
+        self.op(u,self.tmp)
+        u.as_numpy[:] *= 3/5
+        u.axpy(1,self.q2)
+        u.axpy(dt/10, self.tmp)
+
+        return dt
+
+def run(Model, space,
+        dt=None, cfl=None, limiter="default",
+        saveStep=None, subsamp=0):
+    endTime  = Model.endTime
+    name     = Model.name
+    tcount   = 0
+    saveTime = saveStep
+
+    stages = 4
+    operator = femDGOperator(Model, space, limiter=limiter, threading=True, parameters=parameters )
+    # stepper  = Heun(operator, cfl)
+    # stepper  = ExplSSP3(stages,operator,cfl)
+    stepper  = ExplSSP4_10(operator,cfl)
+
+    u_h = space.interpolate(Model.initial, name='u_h')
+    operator.applyLimiter( u_h )
+    vtk = grid.sequencedVTK(name, pointdata=[u_h], subsampling=subsamp)
+    vtk()
+    fixedDt = dt is not None
+
+    while operator.time < endTime:
+        dt = stepper(u_h)
+        # check that solution is meaningful
+        assert not math.isnan( u_h.scalarProductDofs( u_h ) )
+        # increment time and time step counter
+        operator.addToTime(dt)
+        tcount += 1
+        print('[',tcount,']','dt = ', dt, 'time = ',operator.time,
+                'dtEst = ',operator.timeStepEstimate,
+                'elements = ',grid.size(0), flush=True )
+        if saveStep is not None and operator.time > saveTime:
+            vtk()
+            saveTime += saveStep
+    print("number of time steps ", tcount,flush=True)
+    vtk()
+    return u_h, operator.time
+
+Model = problem()
+grid  = structuredGrid(*Model.domain)
+grid.hierarchicalGrid.globalRefine(2)
+space = dgonb(grid, order=3, dimRange=Model.dimRange)
+
+fixedDt = None
+u_h, endTime = run(Model, space, limiter=None, dt=fixedDt, saveStep=0.01)
+
+exact = Model.exact
+if exact is not None:
+    tc = Constant(0, "time")
+    u  = exact(tc)
+    tc.value = endTime # there is an issue here with the generated code in the next step involving the actual endTime
+    error = integrate( grid, dot(u_h-u,u_h-u), order=5 )
+    error = math.sqrt(error)
+    print("error:", error)

pydemo/scalar/scalar.py

 from ufl import *
 from dune.ufl import Space
 import dune.fem
-from dune.femdg.testing import Parameters
-# Parameters are: "Model", "initial", "domain", "endTime", "name", "exact"
 
 class Transport1D:
     dimRange = 1
@@ -38,8 +36,8 @@ def LinearAdvectionDiffusion1D(v,eps):
             # TODO: this method is used in an IMEX method allough the
             # diffusion is implicit - should we change that behaviour?
             # commented out for test of IMEX
-            # def maxDiffusion(t,x,U):
-            #    return eps
+            def maxDiffusion(t,x,U):
+               return eps
         def physical(U):
             return 1
         def jump(U,V):
@@ -101,40 +99,57 @@ def riemanProblem(x,x0,UL,UR):
     return conditional(x<x0,as_vector(UL),as_vector(UR))
 
 def constantTransport():
-    return Parameters(Model=Transport1D, initial=as_vector( [0.1] ),
-                      domain=[[-1, 0], [1, 0.1], [50, 7]], endTime=0.1,
-                      name="constant")
+    Model=Transport1D
+    Model.initial=as_vector( [0.1] )
+    Model.domain=[[-1, 0], [1, 0.1], [50, 7]]
+    Model.endTime=0.1
+    Model.name="constant"
+    return Model
+
 def shockTransport():
-    return Parameters(Model=Transport1D, initial=riemanProblem(x[0],-0.5, [1], [0]),
-                      domain=[[-1, 0], [1, 0.1], [50, 7]], endTime=1.0,
-                      name="shockTransport")
+    Model=Transport1D
+    Model.initial=riemanProblem(x[0],-0.5, [1], [0])
+    Model.domain=[[-1, 0], [1, 0.1], [50, 7]]
+    Model.endTime=1.0,
+    Model.name="shockTransport"
+    return Model
 
 def sinAdvDiffProblem():
     eps = 0.5
     v   = [4,0]
     u0 = lambda t,x: as_vector( [sin(2*pi*(x[0]-t*v[0]))*exp(-t*eps*(2*pi)**2)] )
-    # return Parameters(Model=LinearAdvectionDiffusion1DMixed(v,eps,u0),
-    return Parameters(Model=LinearAdvectionDiffusion1DPeriodic(v,None),
-                      initial=u0(0,x), domain=[[-1, 0], [1, 0.1], [50, 7]], endTime=0.2,
-                      name="sin", exact=lambda t: u0(t,x))
+    Model = Model=LinearAdvectionDiffusion1DMixed(v,eps,u0)
+    Model.initial = u0(0,x)
+    Model.domain = [[-1, 0], [1, 0.1], [50, 7]]
+    Model.endTime = 0.01
+    Model.name = "sin"
+    Model.exact = lambda t: u0(t,x)
+    return Model
 
 def sinDiffProblem():
     eps = 0.5
     u0 = lambda t,x: as_vector( [sin(2*pi*x[0])*exp(-t*eps*(2*pi)**2)] )
-    return Parameters(Model=LinearAdvectionDiffusion1DMixed(None,eps,u0), initial=u0(0,x),
-                      domain=[[-1, 0], [1, 0.1], [50, 7]], endTime=0.2,
-                      name="sin", exact=lambda t: u0(t,x))
+    Model=LinearAdvectionDiffusion1DMixed(None,eps,u0)
+    Model.initial=u0(0,x)
+    Model.domain=[[-1, 0], [1, 0.1], [50, 7]]
+    Model.endTime=0.2
+    Model.name="sin"
+    Model.exact=lambda t: u0(t,x)
+    return Model
 
 def sinTransportProblem():
     v   = [1,0]
     u0 = lambda t,x: as_vector( [sin(2*pi*(x[0]-t*v[0]))] )
-    return Parameters(Model=LinearAdvectionDiffusion1DDirichlet(v,None,u0),
-    # return Parameters(Model=LinearAdvectionDiffusion1DMixed(v,None,u0),
-    # return Parameters(Model=LinearAdvectionDiffusion1D(v,None),
-    # return Parameters(Model=LinearAdvectionDiffusion1DPeriodic(v,None),
-                      initial=u0(0,x),
-                      domain=[[-1, 0], [1, 0.1], [50, 7]], endTime=0.5,
-                      name="sin", exact=lambda t: u0(t,x))
+    Model=LinearAdvectionDiffusion1DDirichlet(v,None,u0)
+    # Model=LinearAdvectionDiffusion1DMixed(v,None,u0)
+    # Model=LinearAdvectionDiffusion1D(v,None)
+    # Model=LinearAdvectionDiffusion1DPeriodic(v,None)
+    Model.initial=u0(0,x)
+    Model.domain=[[-1, 0], [1, 0.2], [25, 7]]
+    Model.endTime=0.5
+    Model.name="sin"
+    Model.exact=lambda t: u0(t,x)
+    return Model
 
 def pulse(eps=None):
     center  = as_vector([ 0.5,0.5 ])
@@ -153,35 +168,46 @@ def pulse(eps=None):
         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 Parameters(Model=LinearAdvectionDiffusion1DDirichlet([ux,uy],eps,u0),
-                      initial=u0(0,x),
-                      domain=[[0, 0], [1, 1], [16,16]], endTime=1.0,
-                      name="pulse"+str(eps), exact=lambda t: u0(t,x))
+    Model=LinearAdvectionDiffusion1DDirichlet([ux,uy],eps,u0)
+    Model.initial=u0(0,x)
+    Model.domain=[[0, 0], [1, 1], [16,16]]
+    Model.endTime=1.0
+    Model.name="pulse"+str(eps)
+    Model.exact=lambda t: u0(t,x)
+    return Model
 def diffusivePulse():
     return pulse(0.001)
 
 def burgersShock():
-    Model = Burgers1D
     UL = 1
     UR = 0
     speed = (UL-UR)/2.
     u0 = lambda t,x: riemanProblem(x[0],-0.5+t*speed,[UL],[UR])
-    return Parameters(Model=Model, initial=u0(0,x),
-                      domain=[[-1, 0], [1, 0.1], [50, 7]], endTime=1.,
-                      name="burgersShock", exact=lambda t: u0(t,x))
+    Model = Burgers1D
+    Model.initial=u0(0,x)
+    Model.domain=[[-1, 0], [1, 0.1], [50, 7]]
+    Model.endTime=1
+    Model.name="burgersShock"
+    Model.exact=lambda t: u0(t,x)
+    return Model
 
 def burgersVW():
     Model = Burgers1D
-    return Parameters(Model=Model, initial=riemanProblem(x[0],0,[-1],[1]),
-                      domain=[[-1, 0], [1, 0.1], [50, 7]], endTime=0.7,
-                      name="burgersShock")
+    Model.initial=riemanProblem(x[0],0,[-1],[1])
+    Model.domain=[[-1, 0], [1, 0.1], [50, 7]]
+    Model.endTime=0.7
+    Model.name="burgersShock"
+    return Model
 
 def burgersStationaryShock():
-    Model = Burgers1D
     u0 = lambda t,x: riemanProblem(x[0],0,[1],[-1])
-    return Parameters(Model=Model, initial=u0(0,x),
-                      domain=[[-1, 0], [1, 0.1], [50, 7]], endTime=0.2,
-                      name="burgersShock", exact=lambda t: u0(t,x))
+    Model = Burgers1D
+    Model.initial=u0(0,x)
+    Model.domain=[[-1, 0], [1, 0.1], [50, 7]]
+    Model.endTime=0.2
+    Model.name="burgersShock"
+    Model.exact=lambda t: u0(t,x)
+    return Model
 
 problems = [ constantTransport, shockTransport,\
              sinDiffProblem, sinTransportProblem, sinAdvDiffProblem,\

python/dune/femdg/_operators.py

@@ -255,10 +255,25 @@ def femDGOperator(Model, space,
 
     _, domainFunctionIncludes, domainFunctionType, _, _, _ = space.storage
     base = 'Dune::Fem::SpaceOperatorInterface< ' + domainFunctionType + '>'
-    return load(includes, typeName,
-                baseClasses = [base],
-                preamble=writer.writer.getvalue()).\
-                Operator( space, advModel, diffModel, parameters=parameters )
+    op = load(includes, typeName,
+             baseClasses = [base],
+             preamble=writer.writer.getvalue()).\
+             Operator( space, advModel, diffModel, parameters=parameters )
+    op._t = t
+    op.time = t.value
+    op.models = [advModel,diffModel]
+    op.space = space
+    def addToTime(self,dt):
+        # print(dt,t.value,self.models[0].time,self.models[1].time)
+        self._t.value += dt
+        self.time += dt
+        self.setTime(self.time)
+        # print(" ",t.value,self.models[0].time,self.models[1].time)
+    def stepTime(self,c,dt):
+        self.setTime(self.time+c*dt)
+    op.addToTime = addToTime.__get__(op)
+    op.stepTime  = stepTime.__get__(op)
+    return op
 
 # RungeKutta solvers
 def rungeKuttaSolver( fullOperator, imex='EX', butchertable=None, parameters={} ):

python/dune/femdg/testing.py

@@ -12,11 +12,20 @@ Parameters = namedtuple("TestingParameters",
                         ["Model", "initial", "domain", "endTime", "name", "exact"])
 Parameters.__new__.__defaults__ = (None,None) # defaults for name, exact
 
-def run(Model, initial, domain, endTime, name, exact,
-        polOrder, limiter="default", startLevel=0,
+def run(Model,
+        initial=None, domain=None, endTime=None, name=None, exact=None, # deprecated - now Model attributes
+        polOrder=1, limiter="default", startLevel=0,
         primitive=None, saveStep=None, subsamp=0,
         dt=None,grid="yasp", space="dgonb", threading=True,
         parameters={}):
+    if initial is not None:
+        print("deprecated usage of run: add initial etc as class atttributes to the Model")
+    else:
+        initial = Model.initial
+        domain = Model.domain
+        endTime = Model.endTime
+        name = Model.name
+        exact = Model.exact
     print("*************************************")
     print("**** Running simulation",name)
     print("*************************************")
@@ -31,10 +40,10 @@ def run(Model, initial, domain, endTime, name, exact,
                 if 2*i+1 in bnd:
                     assert(2*i+2 in bnd)
                     periodic[i] = False
-        print("Setting periodic boundaries",periodic,flush=True)
         # create domain and grid
         domain   = cartesianDomain(x0,x1,N,periodic=periodic,overlap=0)
         grid     = create.grid(grid,domain)
+        print("Setting periodic boundaries",periodic,flush=True)
     except TypeError: # assume the 'domain' is already a gridview
         grid = domain
     # initial refinement of grid