From 84e2bda61eb18dd436fd6125c99df3ebf2312ef7 Mon Sep 17 00:00:00 2001
From: dedner <a.s.dedner@warwick.ac.uk>
Date: Tue, 31 Jul 2018 20:27:46 +0100
Subject: [PATCH] set up a heat equation example (not working)

---
 pydemo/scalar/scalar.py     | 17 ++++++++++-------
 pydemo/scalar/shock_tube.py | 16 +++++++++-------
 2 files changed, 19 insertions(+), 14 deletions(-)

diff --git a/pydemo/scalar/scalar.py b/pydemo/scalar/scalar.py
index e322002d..4c92d81b 100644
--- a/pydemo/scalar/scalar.py
+++ b/pydemo/scalar/scalar.py
@@ -26,7 +26,7 @@ class LinearAdvectionDiffusion1D:
     # def F_c(t,x,U):
     #     return as_matrix( [ [U[0], 0] ] )
     def F_v(t,x,U,DU):
-        return 0.1*DU
+        return 0.5*DU
     def maxLambda(t,x,U,n):
         return abs(n[0])
     def velocity(t,x,U):
@@ -35,11 +35,13 @@ class LinearAdvectionDiffusion1D:
         return 1
     def jump(U,V):
         return U-V
-class LinearAdvectionDiffusion1DDirichlet(LinearAdvectionDiffusion1D):
-    def dirichletValue(t,x,u):
-        return as_vector( [cos(2*pi*x[0])] )
-    boundaryValue = {}
-    for i in range(1,5): boundaryValue.update( {i: dirichletValue} )
+def LinearAdvectionDiffusion1DDirichlet(bnd):
+    class Model(LinearAdvectionDiffusion1D):
+        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
@@ -77,7 +79,8 @@ def shockTransport():
     return Transport1D, riemanProblem(x[0],-0.5, [1], [0])
 
 def sinProblem():
-    return LinearAdvectionDiffusion1DDirichlet, as_vector( [sin(2*pi*x[0])] )
+    u0 = as_vector( [sin(2*pi*x[0])] )
+    return LinearAdvectionDiffusion1DDirichlet(u0), u0
 
 def cosProblem():
     return LinearAdvectionDiffusion1DNeuman, as_vector( [cos(2*pi*x[0])] )
diff --git a/pydemo/scalar/shock_tube.py b/pydemo/scalar/shock_tube.py
index 0f66866d..91e2a795 100644
--- a/pydemo/scalar/shock_tube.py
+++ b/pydemo/scalar/shock_tube.py
@@ -5,9 +5,9 @@ import dune.create as create
 from dune.models.elliptic.formfiles import loadModels
 from dune.femdg import createFemDGSolver
 
-from scalar import shockTransport as problem
+# from scalar import shockTransport as problem
 # from scalar import constantTransport as probelm
-# from scalar import cosProblem as problem
+from scalar import sinProblem as problem
 # from scalar import cosProblem as problem
 # from scalar import burgersShock as problem
 # from scalar import burgersVW as problem
@@ -15,33 +15,35 @@ from scalar import shockTransport as problem
 
 Model, initial = problem()
 
-
+parameter.append({"fem.verboserank": 0})
 parameter.append("parameter")
-parameter.append({"fem.verboserank": -1})
 
 x0,x1,N = [-1, 0], [1, 0.1], [50, 7]
 grid = structuredGrid(x0,x1,N)
+
 # grid = create.grid("ALUSimplex", cartesianDomain(x0,x1,N))
 dimR     = 1
 t        = 0
+dt       = 0.001
 saveStep = 0.01
 saveTime = saveStep
 count    = 0
 endTime  = 1.
 
 def useODESolver():
-    global count, t, saveTime
+    global count, t, dt, saveTime
 
     spaceName = "dgonb"
     polOrder = 2
     space = create.space(spaceName, grid, order=polOrder, dimrange=dimR)
     u_h   = space.interpolate(initial, name='u_h')
-    operator = createFemDGSolver( Model, space )
+    operator = createFemDGSolver( Model, space, limiter=None )
+
+    operator.setTimeStepSize(dt)
 
     start = time.time()
     grid.writeVTK(Model.name, celldata=[u_h], number=count)
     while t < endTime:
-        # operator.applyLimiter( u_h );
         operator.solve(target=u_h)
         dt = operator.deltaT()
         t += dt
-- 
GitLab