Merge branch 'feature/scaling-limiter' into feature/dg-operator

pydemo/euler/euler.ufl

+gamma = 1.4
+class Compressible2DEuler:
+    def pressure(U):
+        return (U[3]-(U[1]**2+U[2]**2)/2)*(gamma-1)
+    def toCons(U):
+        return [U[0],U[0]*U[1],U[0]*U[2],U[3]/(gamma-1)+U[0]/2*(U[1]**2+U[2]**2)]
+    def toPrim(U):
+        return [U[0],U[1]/U[0],U[2]/U[0],Compressible2DEuler.pressure(U)]
+    def F_c(U):
+        rho, u,v, p = Compressible2DEuler.toPrim(U)
+        rE = U[-1]
+        res = as_matrix( [ [rho*u, rho*v],
+                           [rho*u**2 + p, rho*u*v],
+                           [rho*u*v, rho*v**2 + p],
+                           [(rE+p)*u, (rE+p)*v]
+                         ])
+        return res
+    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
+
+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)

pydemo/euler/shock_tube.py

+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
+
+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)
+
+import euler
+
+Euler = euler.Model
+NumFLux = euler.NumFlux
+
+u = TrialFunction(space)
+v = TestFunction(space)
+dt = NamedConstant(triangle, "dt")    # time step
+x = SpatialCoordinate(space.cell())
+n = FacetNormal(space.cell())
+
+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)
+
+u_h = space.interpolate(initial, name='u_h')
+u_h_n = u_h.copy(name="previous")
+
+# 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.model.dt = 5e-4
+
+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)