minor change.

pydemo/euler/paper.py

@@ -18,11 +18,12 @@ class Model:
     # interface methods for model
     def F_c(t,x,U):
         rho, v, p = Model.toPrim(U)
-        return as_matrix( [
-                  [rho*v[0], rho*v[1]],
-                  [rho*v[0]*v[0] + p, rho*v[0]*v[1]],
-                  [rho*v[0]*v[1], rho*v[1]*v[1] + p],
-                  [(U[3]+p)*v[0], (U[3]+p)*v[1]] ] )
+        rE = U[3]
+        v = numpy.array(v)
+        res = numpy.vstack([ rho*v,
+                             rho*numpy.outer(v,v) + p*numpy.eye(2),
+                             (rE+p)*v ])
+        return as_matrix(res)
     # simple 'outflow' boundary conditions on all boundaries
     boundary = {range(1,5): lambda t,x,U: U}
 
@@ -33,22 +34,21 @@ class Model:
 
 ########################################################################
 
-print("\n Part 1\n")
-
-# Part 1: basic set up and time loop - no limiter and fixed time step
-#         using constant initial data
-gridView = structuredGrid([-1,-1],[1,1],[40,40])
-space = dgonb( gridView, order=3, dimRange=4)
-operator = femDGOperator(Model, space, limiter=None)
-stepper  = femdgStepper(order=3, operator=operator)
-u_h = space.interpolate([1.4,0,0,1], name='u_h')
-t  = 0
-# don't show vtk in paper
-vtk = gridView.sequencedVTK("paperA", pointdata=[u_h])
-while t < 0.1:
-    operator.setTime(t)
-    t += stepper(u_h, dt=0.001)
-vtk()
+if False:
+    # Part 1: basic set up and time loop - no limiter and fixed time step
+    #         using constant initial data
+    gridView = structuredGrid([-1,-1],[1,1],[40,40])
+    space = dgonb( gridView, order=3, dimRange=4)
+    operator = femDGOperator(Model, space, limiter=None)
+    stepper  = femdgStepper(order=3, operator=operator)
+    u_h = space.interpolate([1.4,0,0,1], name='u_h')
+    t  = 0
+    # don't show vtk in paper
+    vtk = gridView.sequencedVTK("paperA", pointdata=[u_h])
+    while t < 0.1:
+        operator.setTime(t)
+        t += stepper(u_h, dt=0.001)
+    vtk()
 
 # Part 2: add limiter and methods for troubled cell indicator
 def velocity(t,x,U):
@@ -66,59 +66,50 @@ Model.velocity = velocity
 Model.physical = physical
 Model.jump     = jump
 
-import os
-fluxHeader = os.getcwd() + "/llf.hh"
+if False:
 
-operator = femDGOperator(Model, space, advectionFlux=fluxHeader, limiter="MinMod")
-stepper  = femdgStepper(order=3, operator=operator)
-x = SpatialCoordinate(space)
-u_h.interpolate( conditional(dot(x,x)<0.1,as_vector([1,0,0,2.5]),
-                                          as_vector([0.125,0,0,0.25])) )
-operator.applyLimiter(u_h)
-t  = 0
-vtk = gridView.sequencedVTK("paperB", pointdata=[u_h])
+    import os
+    fluxHeader = os.getcwd() + "/llf.hh"
 
-print("Start Part 2\n")
-while t < 0.4:
+    operator = femDGOperator(Model, space, advectionFlux=fluxHeader, limiter="MinMod")
+    stepper  = femdgStepper(order=3, operator=operator)
+    x = SpatialCoordinate(space)
+    u_h.interpolate( conditional(dot(x,x)<0.1,as_vector([1,0,0,2.5]),
+                                              as_vector([0.125,0,0,0.25])) )
+    operator.applyLimiter(u_h)
+    t  = 0
+    vtk = gridView.sequencedVTK("paperB", pointdata=[u_h])
+    while t < 0.4:
+        vtk()
+        operator.setTime(t)
+        t += stepper(u_h)
     vtk()
-    operator.setTime(t)
-    t += stepper(u_h)
-vtk()
-
-# Part 3: FV on polygonal grid
-from dune.generator import algorithm
-from dune.polygongrid import voronoiDomain, polygonGrid
-from dune.grid import reader
-
-import numpy
-
-x0 = [-1,-1]
-x1 = [ 1, 1]
-N  = [40,40]
-#boundingBox = numpy.array([ x0, x1 ])
-#vb = voronoiDomain(N[0]*N[1], boundingBox, seed=1234)
 
-domain = (reader.dgf, "triangle.dgf")
-gridView = polygonGrid( domain, dualGrid=False )
-fvspc = finiteVolume( gridView, dimRange=4)
-space = dgonb( gridView, order=1, dimRange=4, caching=False )
-u_h   = space.interpolate( conditional(dot(x,x)<0.1,as_vector([1,0,0,2.5]),
-                                                    as_vector([0.125,0,0,0.25])),
-                                                    name="uh")
-fvU = fvspc.interpolate( [0,0,0,0], name = "fvU" )
+if True:
+    # Part 3: FV on polygonal grid
+    from dune.generator import algorithm
+    from dune.polygongrid import voronoiDomain, polygonGrid
+    import numpy
 
-operator = femDGOperator(Model, space, limiter="MinMod")
-stepper  = femdgStepper(order=1, operator=operator, cfl=0.4)
-operator.applyLimiter(u_h)
-count = 0
-gridView.writeVTK("paperC", celldata=[fvU], number=count)
-t  = 0
-print("Start Part 3\n")
-while t < 0.4:
-    assert not math.isnan( u_h.scalarProductDofs( u_h ) )
-    operator.setTime(t)
-    t += stepper(u_h)
-    print(t)
-    fvU.interpolate( u_h )
-    count += 1
-    gridView.writeVTK("paperC", celldata=[fvU], number=count)
+    x0 = [-1,-1]
+    x1 = [ 1, 1]
+    N  = [40,40]
+    boundingBox = numpy.array([ x0, x1 ])
+    vb = voronoiDomain(N[0]*N[1], boundingBox, seed=1234)
+    gridView = polygonGrid( vb )
+    # space = finiteVolume( gridView, dimRange=4)
+    space = dgonb( gridView, order=2, dimRange=4)
+    x = SpatialCoordinate(space)
+    u_h   = space.interpolate( conditional(dot(x,x)<0.1,as_vector([1,0,0,2.5]),
+                                                        as_vector([0.125,0,0,0.25])),
+                                                        name="uh")
+    operator = femDGOperator(Model, space, limiter="MinMod")
+    stepper  = femdgStepper(order=3, operator=operator, cfl=0.4)
+    operator.applyLimiter(u_h)
+    t  = 0
+    while t < 0.4:
+        assert not math.isnan( u_h.scalarProductDofs( u_h ) )
+        operator.setTime(t)
+        t += stepper(u_h)
+        print(t)
+    algorithm.run('vtkout', 'vtkout.hh', u_h )