growing domain turing patterns with extend

growing_domain/growing_domain_patterns/growing_domain_patterns.py

@@ -20,9 +20,9 @@ class GrowingDomainPatterns(object):
         #growth_direction = ('x[0]', )
         self.growth_direction = df.Expression(growth_direction, degree=1)
         self.sigma = df.Expression(sigma[self.growth], 
-                                        L0=self.system_size,
-                                        b=self.growth_parameter, 
-                                        t=0, degree=1)
+                                        L0 = self.system_size,
+                                        b = self.growth_parameter, 
+                                        t = 0, degree=1)
         # set up mesh
         if self.dimension==1:
             self.mesh = df.IntervalMesh(self.resolution, 0, self.system_size)
@@ -35,6 +35,7 @@ class GrowingDomainPatterns(object):
         scalar_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1)
         mixed_element = df.MixedElement([scalar_element, scalar_element])
         self.function_space = df.FunctionSpace(self.mesh, mixed_element)
+
         self.f = df.Function(self.function_space)
         self.c1, self.c2 = df.split(self.f)
         self.f0 = df.Function(self.function_space)
@@ -44,7 +45,6 @@ class GrowingDomainPatterns(object):
         self.VFS = df.VectorFunctionSpace(self.mesh, 'P', 1)        
         self.velocity = df.project(self.sigma * self.growth_direction, self.VFS)
 
-
         c1form        = (df.inner((self.c1 - self.c10)/self.timestep, tc1) 
                         + df.inner(df.div(self.velocity*self.c10), tc1) 
                         + self.Dc1 * df.inner(df.nabla_grad(self.c1), df.nabla_grad(tc1))
@@ -57,14 +57,21 @@ class GrowingDomainPatterns(object):
 
         self.form = (c1form + c2form)*df.dx(self.mesh)
 
-
-    def solve(self):
+    def solve(self, extend):
         savesteps = int(self.savetime/self.timestep)           
         maxsteps = int(self.maxtime/self.timestep)
         self.time = 0.0
-        c1File = df.XDMFFile('%s_c1.xdmf' %params['timestamp'])   # create the file here
-        c2File = df.XDMFFile('%s_c2.xdmf' %params['timestamp']) 
-
+        self.c1File = df.XDMFFile('%s_c1.xdmf' %params['timestamp'])   # create the file here
+        self.c2File = df.XDMFFile('%s_c2.xdmf' %params['timestamp']) 
+
+        if extend:
+            SFS = df.FunctionSpace(self.mesh,'P',1)
+            c10, c20 = df.Function(SFS), df.Function(SFS)
+            self.c1File.read_checkpoint(c10, 'c1', -1)
+            self.c2File.read_checkpoint(c20, 'c2', -1)
+            c10 = df.interpolate(c10, self.function_space.sub(0).collapse())
+            c20 = df.interpolate(c20, self.function_space.sub(1).collapse())
+        else:
             # initial condition
             c10 = df.interpolate(df.Constant(self.a), self.function_space.sub(0).collapse())
             noise_c1 = self.noise_level * (2*np.random.random(c10.vector().size())-1)
@@ -77,8 +84,8 @@ class GrowingDomainPatterns(object):
         df.assign(self.f0, [c10, c20])
 
         # save data
-        c1File.write_checkpoint(c10, 'c1', self.time, append=True)
-        c2File.write_checkpoint(c20, 'c2', self.time, append=True)
+        self.c1File.write_checkpoint(c10, 'c1', self.time, append=True)
+        self.c2File.write_checkpoint(c20, 'c2', self.time, append=True)
 
         # time stepping
         for steps in progressbar.progressbar(range(1, maxsteps + 1)):
@@ -93,36 +100,47 @@ class GrowingDomainPatterns(object):
             c1, c2 = self.f0.split(deepcopy=True)
 
             if steps % savesteps == 0:
-                c1File.write_checkpoint(c1, 'c1', self.time, append=True)
-                c2File.write_checkpoint(c2, 'c2', self.time, append=True)
+                self.c1File.write_checkpoint(c1, 'c1', self.time, append=True)
+                self.c2File.write_checkpoint(c2, 'c2', self.time, append=True)
 
             # move mesh
             displacement = df.project(self.velocity*self.timestep, self.VFS)
             df.ALE.move(self.mesh, displacement)
             self.time += self.timestep
 
-        c1File.close()
-        c2File.close()
+        self.c1File.close()
+        self.c2File.close()
     
 
 if __name__ == '__main__':
     import json, datetime
-    assert os.path.isfile('parameters.json'), 'parameters.json file not found'  # assert that parameters.json is a valid file, otherwise 
+    import argparse
                                                                                     # give an error message parameters.json file not found
+    parser = argparse.ArgumentParser()
+    parser.add_argument('-j','--jsonfile', help='json file', default='growing_domain_patterns_parameters.json')
+    parser.add_argument('-t','--time', help='time to extend', type=float, default=50)
+    args = parser.parse_args()
 
-    # load the parameters
-    with open('parameters.json') as jsonFile:                    
+    assert os.path.isfile(args.jsonfile), '%s file not found' % args.jsonfile
+    with open(args.jsonfile) as jsonFile:
         params = json.load(jsonFile)
 
-    # parse parameters
-    assert params['dimension'] in (1,2)
-    assert params['growth'] in ('none','linear','exponential'), 'Unknown growth model'
-    
+    if args.jsonfile=='growing_domain_patterns_parameters.json':
+        print('Fresh run')
+        extend = False
         timestamp = datetime.datetime.now().strftime("%d%m%y-%H%M%S")
+        # current timestamp is the jobID
         params['timestamp'] = timestamp
 
-    gdp = GrowingDomainPatterns(params)
-    gdp.solve()
+    else:
+        print('Extending %s' % params['timestamp'])
+        extend = True
+        oldMaxTime = params['maxtime']
+        params['maxtime'] = args.time
         
+    gdp = GrowingDomainPatterns(params)
+    gdp.solve(extend)
+    if extend:
+        params['maxtime'] = oldMaxTime + params['maxtime']
     with open(params['timestamp'] + '_parameters.json', "w") as fp:
         json.dump(params, fp, indent=4)    
\ No newline at end of file