testing save and extend on a simple fixed domain diffusion problem

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
+
+class diffusion(object):
+    def __init__(self, parameters):
+        # read in parameters
+        for key in parameters:
+            setattr(self, key, parameters[key])
+        
+        self.mesh = df.IntervalMesh(self.resolution, 
+                                    -self.L/2, 
+                                     self.L/2)
+
+        self.function_space = df.FunctionSpace(self.mesh,'P', 1)
+        self.c = df.Function(self.function_space)
+        self.tc = df.TestFunction(self.function_space)
+        self.c0 = df.Function(self.function_space)
+        
+    def setup_initial_conditions(self, initc):
+        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])
+
+    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):
+        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)
+        self.setup_weak_forms()
+
+        # if fresh run, write initial condition
+        if not extend:
+            self.cFile.write_checkpoint(self.c0, 'concentration', self.time)
+
+        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.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,
+    "L": 1.0, 
+    "D": 0.1,
+    "timestep": 0.1,
+    "maxtime": 0.5
+
+}   
\ 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',
+        default='parameters.json')
+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
+
+# load the parameters
+with open(args.jsonfile) as jsonFile:                    
+    parameters = json.load(jsonFile)
+
+if args.jsonfile=='parameters.json':
+    extend = False
+    print('Fresh run...')
+    timestamp = datetime.datetime.now().strftime("%d%m%y-%H%M%S")
+    parameters['timestamp'] = timestamp
+    parametersFile = timestamp + '_parameters.json'
+    initc = None
+    initTime = 0.0
+
+else:
+    extend = True
+    print('Extending run')
+    parametersFile = args.jsonfile
+
+    steps = int(parameters['maxtime']/parameters['timestep'])
+    
+    mesh = df.IntervalMesh(parameters['resolution'], 
+                            - 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)
+    cFile   = df.XDMFFile('%s_concentration.xdmf' % parameters['timestamp'])
+    cFile.read_checkpoint(initc, 'concentration', steps)
+    cFile.close()
+    initTime = parameters['maxtime']
+    parameters['maxtime'] = args.time
+
+diff = diffusion(parameters)
+diff.solve(initc, initTime, extend)
+
+if extend:
+    parameters['maxtime'] = initTime + parameters['maxtime']
+
+with open(parametersFile, "w") as jsonFile:
+    json.dump(parameters, jsonFile, indent=4, sort_keys=True)
+
+from viz_extend_for_diffusion import visualize
+visualize(parameters, DIR='')
\ No newline at end of file

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'])
+    times = np.arange(savesteps+1) * params['timestep']
+    
+    mesh = df.IntervalMesh(params['resolution'], 
+                                - params['L']/2, 
+                                  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']))
+
+    # Reading data
+    print('Reading data...')
+    for steps in progressbar.progressbar(range(savesteps+1)):
+        SFS = df.FunctionSpace(mesh, 'P', 1)
+        ci= df.Function(SFS)
+            
+        cFile.read_checkpoint(ci, 'concentration', steps)
+        c[steps] = ci.compute_vertex_values(mesh)
+
+    cFile.close()
+
+    # interactive plot
+    plt.rcParams.update({'font.size': 15})
+
+    fig, axes = plt.subplots(1, 1, sharex=True, figsize=(8,8))
+    axes.set_xlabel(r'$x$')
+    axes.set_ylabel(r'$c$')
+    x = mesh.coordinates()[:,0]
+    cplot, = axes.plot(x, c[0], 'g-', ms=3)
+    
+    def update(value):
+        ti = np.abs(times-value).argmin()
+        cplot.set_ydata(c[ti])
+        plt.draw()
+        
+    sax = plt.axes([0.1, 0.92, 0.7, 0.02])
+    slider = Slider(sax, r'$t/\tau$', min(times), max(times),
+            valinit=min(times), valfmt='%3.1f',
+            fc='#999999')
+    slider.drawon = False
+    slider.on_changed(update)
+    plt.show()
+
+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:
+        params = json.load(jsonFile)
+
+    visualize(params, DIR=os.path.dirname(args.jsonfile))
\ No newline at end of file