From e9ba2c0c3020d781f1a80ddd51e0565e97507289 Mon Sep 17 00:00:00 2001
From: dedner <a.s.dedner@warwick.ac.uk>
Date: Fri, 3 Aug 2018 11:52:28 +0100
Subject: [PATCH] added error computation for pulse problem - for other trst
 cases exact solution needs to be added

---
 pydemo/scalar/paramSolver   |  6 +++---
 pydemo/scalar/scalar.py     |  2 +-
 pydemo/scalar/shock_tube.py | 14 ++++++++------
 3 files changed, 12 insertions(+), 10 deletions(-)

diff --git a/pydemo/scalar/paramSolver b/pydemo/scalar/paramSolver
index 08adec49..a2df1b67 100644
--- a/pydemo/scalar/paramSolver
+++ b/pydemo/scalar/paramSolver
@@ -15,9 +15,9 @@ fem.timeprovider.updatestep: 1
 # parameter for the implicit solvers
 # fem.differenceoperator.eps: 1e-12
 fem.solver.gmres.restart: 50
-fem.solver.newton.verbose: true
-fem.solver.newton.linear.verbose: true
-fem.solver.verbose: true
+fem.solver.newton.verbose: false
+fem.solver.newton.linear.verbose: false
+fem.solver.verbose: false
 fem.solver.newton.maxlineariterations: 100000
 fem.solver.newton.tolerance: 1e-07
 fem.solver.newton.linabstol: 1.25e-9
diff --git a/pydemo/scalar/scalar.py b/pydemo/scalar/scalar.py
index 446b25f2..63580c72 100644
--- a/pydemo/scalar/scalar.py
+++ b/pydemo/scalar/scalar.py
@@ -117,7 +117,7 @@ def pulse():
         yq = (-x0*sin(4.0*t) + x1*cos(4.0*t))
         return as_vector( [(sig2/ (sig2PlusDt4) ) * exp (-( xq*xq + yq*yq ) / sig2PlusDt4 )] )
     return LinearAdvectionDiffusion1DDirichlet([ux,uy],eps,u0), u0(0,x),\
-           [0, 0], [1, 1], [8,8], 1.0,\
+           [0, 0], [1, 1], [16,16], 1.0,\
            "pulse", lambda t: u0(t,x)
 
 def burgersShock():
diff --git a/pydemo/scalar/shock_tube.py b/pydemo/scalar/shock_tube.py
index 591027b1..2ee9ad80 100644
--- a/pydemo/scalar/shock_tube.py
+++ b/pydemo/scalar/shock_tube.py
@@ -1,4 +1,4 @@
-import time
+import time, math
 from dune.grid import structuredGrid, cartesianDomain
 from dune.fem import parameter
 from dune.fem.function import integrate
@@ -9,8 +9,8 @@ from ufl import dot
 
 # from scalar import shockTransport as problem
 # from scalar import constantTransport as probelm
-from scalar import sinProblem as problem
-# from scalar import pulse as problem
+# from scalar import sinProblem as problem
+from scalar import pulse as problem
 
 Model, initial, x0,x1,N, endTime, name, exact = problem()
 
@@ -22,8 +22,8 @@ grid = structuredGrid(x0,x1,N)
 # grid = create.grid("ALUSimplex", cartesianDomain(x0,x1,N))
 dimR     = 1
 t        = 0
-dt       = 0.001
-saveStep = 0.01
+dt       = 0.005
+saveStep = 0.05
 saveTime = saveStep
 count    = 0
 
@@ -51,6 +51,8 @@ def useODESolver():
             saveTime += saveStep
     grid.writeVTK(name, celldata=[u_h], number=count, subsampling=2)
     print("time loop:",time.time()-start)
-    # error = integrate( dot(u_h-exact(t),u_h-exact(t)), order=5 )
+    if exact is not None:
+        error = integrate( grid, dot(u_h-exact(t),u_h-exact(t)), order=5 )
+        print("error:", math.sqrt(error) )
 
 useODESolver()
-- 
GitLab