project the IC, not interpolate. Computing errornorm

simple_extend_code/111223-235806_parameters.json

+{
+    "D": 0.1,
+    "L": 1.0,
+    "maxtime": 1.0,
+    "resolution": 128,
+    "timestamp": "111223-235806",
+    "timestep": 0.1
+}
\ No newline at end of file

simple_extend_code/111223-235856_parameters.json

+{
+    "D": 0.1,
+    "L": 1.0,
+    "maxtime": 1,
+    "resolution": 128,
+    "timestamp": "111223-235856",
+    "timestep": 0.1
+}
\ No newline at end of file

simple_extend_code/121223-065829_parameters.json

+{
+    "D": 0.1,
+    "L": 1.0,
+    "maxtime": 1.0,
+    "resolution": 128,
+    "timestamp": "121223-065829",
+    "timestep": 0.1
+}
\ No newline at end of file

simple_extend_code/121223-065916_parameters.json

+{
+    "D": 0.1,
+    "L": 1.0,
+    "maxtime": 1,
+    "resolution": 128,
+    "timestamp": "121223-065916",
+    "timestep": 0.1
+}
\ No newline at end of file

simple_extend_code/compare.py

+import dolfin as df
+import numpy as np
+
+c1File = df.XDMFFile('121223-065916_concentration.xdmf')
+c2File = df.XDMFFile('121223-065829_concentration.xdmf')
+
+mesh = df.IntervalMesh(128, -0.5, 0.5)
+times = np.arange(0.0, 1.1, 0.1)
+SFS = df.FunctionSpace(mesh, 'P', 1)
+
+for steps in range(len(times)):
+    c1 = df.Function(SFS)
+    c2 = df.Function(SFS)
+    c1File.read_checkpoint(c1, 'concentration', steps)
+    c2File.read_checkpoint(c2, 'concentration', steps)
+    print(df.assemble((c2 - c1)**2 * df.dx(mesh))**0.5)
+
+c1File.close()
+c2File.close()

simple_extend_code/extend_for_diffusion.py

 import numpy as np
 import dolfin as df
 import progressbar
-import os
 
 df.set_log_level(df.LogLevel.ERROR)
 df.parameters['form_compiler']['optimize'] = True
@@ -25,21 +24,19 @@ class diffusion(object):
         if initc is None:
             self.c0.interpolate(df.Expression('1 + 0.2 * cos(pi*x[0]/L)', pi=np.pi, L=self.L, degree=1))
         else:
-            self.c0.interpolate(initc)
-        print('c0[-1]=', self.c0.compute_vertex_values(self.mesh)[-1])
+            self.c0 = df.project(initc, self.function_space)
 
     def setup_weak_forms(self):
         self.form = (df.inner((self.c - self.c0)/self.timestep, self.tc) 
             + self.D * df.inner(df.nabla_grad(self.c), df.nabla_grad(self.tc))
             ) * df.dx
 
-
     def solve(self, initc=None, initTime=0.0, extend=False):
+        # create file if it doesn't exist, open if it exists
         self.cFile = df.XDMFFile('%s_concentration.xdmf' % self.timestamp)
         
         # time-variables
         self.time = initTime
-        print('time=', self.time)
         maxsteps = int(self.maxtime/self.timestep)
 
         self.setup_initial_conditions(initc)
@@ -51,12 +48,9 @@ class diffusion(object):
 
         for steps in progressbar.progressbar(range(1, maxsteps + 1)):
             self.time += self.timestep        
-            print('time=', self.time)    
             df.solve(self.form == 0, self.c)
-            print('c[-1]=', self.c.compute_vertex_values(self.mesh)[-1])
+            self.cFile.write_checkpoint(self.c, 'concentration', self.time, append=True)
             self.c0.assign(self.c)
-            self.cFile.write_checkpoint(self.c, 'concentration', 
-                                            self.time, append=True)
 
         self.cFile.close()
   

simple_extend_code/parameters.json

 {
-    "resolution": 1000,
+    "resolution": 128,
     "L": 1.0, 
     "D": 0.1,
     "timestep": 0.1,
-    "maxtime": 0.5
+    "maxtime": 1
+
 
 }   
\ No newline at end of file

simple_extend_code/run_extend_for_diffusion.py

 import json, datetime
 import os
 from extend_for_diffusion import diffusion
-from viz_extend_for_diffusion import visualize
 import argparse
 import dolfin as df
-import h5py
 
 parser = argparse.ArgumentParser()
-parser.add_argument('-j','--jsonfile', help='data file',
+parser.add_argument('-j','--jsonfile', help ='data file',
         default='parameters.json')
-parser.add_argument('-t','--time', help='time to run', type=float)
+parser.add_argument('-t','--time', help ='time to run', type=float)
 args = parser.parse_args()
 
 assert os.path.isfile(args.jsonfile), '%s file not found' % args.jsonfile
@@ -18,7 +16,7 @@ assert os.path.isfile(args.jsonfile), '%s file not found' % args.jsonfile
 with open(args.jsonfile) as jsonFile:                    
     parameters = json.load(jsonFile)
 
-if args.jsonfile=='parameters.json':
+if args.jsonfile == 'parameters.json':
     extend = False
     print('Fresh run...')
     timestamp = datetime.datetime.now().strftime("%d%m%y-%H%M%S")
@@ -29,7 +27,7 @@ if args.jsonfile=='parameters.json':
 
 else:
     extend = True
-    print('Extending run')
+    print('Extending run...')
     parametersFile = args.jsonfile
 
     steps = int(parameters['maxtime']/parameters['timestep'])
@@ -38,10 +36,6 @@ else:
                             - parameters['L']/2, 
                               parameters['L']/2)
     
-    # read geometry from h5 file
-    var = 'concentration'
-    h5 = h5py.File( '%s_%s.h5' % (parameters['timestamp'], var,), 'r+')
-
     # Read data
     SFS = df.FunctionSpace(mesh, 'P', 1)
     initc = df.Function(SFS)

simple_extend_code/viz_extend_for_diffusion.py

 import dolfin as df
 import numpy as np
 import os
-import h5py
 from matplotlib.widgets import Slider
 import progressbar
 import matplotlib.pyplot as plt
-from tempfile import TemporaryDirectory
 
 def visualize(params, DIR=''):
     savesteps = int(params['maxtime']/params['timestep'])
@@ -16,16 +14,13 @@ def visualize(params, DIR=''):
                                   params['L']/2)
     # Read data
     c = np.zeros((len(times), mesh.num_vertices()))
-
-    cFile   = df.XDMFFile(os.path.join(DIR, 
-                                    '%s_concentration.xdmf' % params['timestamp']))
+    cFile   = df.XDMFFile(os.path.join(DIR, '%s_concentration.xdmf' % params['timestamp']))
 
     # Reading data
     print('Reading data...')
     for steps in progressbar.progressbar(range(savesteps+1)):
         SFS = df.FunctionSpace(mesh, 'P', 1)
-        ci= df.Function(SFS)
-            
+        ci = df.Function(SFS)            
         cFile.read_checkpoint(ci, 'concentration', steps)
         c[steps] = ci.compute_vertex_values(mesh)
 
@@ -57,7 +52,6 @@ if __name__ == '__main__':
     import argparse, json
     parser = argparse.ArgumentParser()
     parser.add_argument('-j','--jsonfile', help='parameters file', required=True)
-    parser.add_argument('--onscreen', action=argparse.BooleanOptionalAction)
     args = parser.parse_args()
 
     with open(args.jsonfile) as jsonFile: