turing patterns on a growing domain

growing_domain/growing_domain_patterns/growing_domain_patterns.py

+import dolfin as df
+import numpy as np
+import progressbar
+import scipy as sc
+import os
+df.set_log_level(df.LogLevel.ERROR)
+df.parameters['form_compiler']['optimize'] = True
+
+class GrowingDomainPatterns(object):
+    def __init__(self, parameters):
+        # read in parameters
+        for key in parameters:
+            setattr(self, key, parameters[key])
+        sigma = {'none': '0.0',
+                'linear': 'b/(L0 + b*t)',
+                'exponential': 'b'}
+
+        # define the growth velocity field
+        growth_direction = tuple(['x['+str(i)+']' for i in range(self.dimension)])  # r hat = (x hat, y hat)
+        #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)
+        # set up mesh
+        if self.dimension==1:
+            self.mesh = df.IntervalMesh(self.resolution, 0, self.system_size)
+        elif self.dimension==2: 
+            self.mesh = df.RectangleMesh(df.Point(0.0,0.0), df.Point(self.system_size, self.system_size), self.resolution, self.resolution)
+        elif self.dimension==3:
+            self.mesh = df.Mesh('sphere.xml.gz')
+        
+        # create function space, define function, test function
+        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)
+        self.c10, self.c20 = df.split(self.f0)
+        tc1, tc2 = df.TestFunctions(self.function_space)
+
+        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))
+                        - df.inner(self.a + self.c1**2*self.c2 - self.b*self.c1 -self.c1, tc1) )
+
+        c2form        = (df.inner((self.c2 - self.c20)/self.timestep, tc2) 
+                        + df.inner(df.div(self.velocity*self.c20), tc2) 
+                        + self.Dc2 * df.inner(df.nabla_grad(self.c2), df.nabla_grad(tc2))
+                        - df.inner(- self.c1**2*self.c2 + self.b*self.c1, tc2) )
+
+        self.form = (c1form + c2form)*df.dx(self.mesh)
+
+
+    def solve(self):
+        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']) 
+
+        # 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)
+        c10.vector()[:] = c10.vector()[:] + noise_c1
+
+        c20 = df.interpolate(df.Constant(self.b/self.a), self.function_space.sub(1).collapse())
+        noise_c2 = self.noise_level * (2*np.random.random(c20.vector().size())-1)
+        c20.vector()[:] = c20.vector()[:] + noise_c2
+
+        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)
+
+        # time stepping
+        for steps in progressbar.progressbar(range(1, maxsteps + 1)):
+            # get velocity
+            self.sigma.t = self.time
+            self.velocity.assign(df.project(self.sigma*self.growth_direction, self.VFS))
+            # solve
+            df.solve(self.form == 0, self.f)
+            # update
+            self.f0.assign(self.f)
+            # save data
+            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)
+
+            # 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()
+    
+
+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 
+                                                                                # give an error message parameters.json file not found
+    
+    # load the parameters
+    with open('parameters.json') as jsonFile:                    
+        params = json.load(jsonFile)
+        
+    # parse parameters
+    assert params['dimension'] in (1,2)
+    assert params['growth'] in ('none','linear','exponential'), 'Unknown growth model'
+    
+    timestamp = datetime.datetime.now().strftime("%d%m%y-%H%M%S")
+    params['timestamp'] = timestamp
+        
+    gdp = GrowingDomainPatterns(params)
+    gdp.solve()
+    
+    with open(params['timestamp'] + '_parameters.json', "w") as fp:
+        json.dump(params, fp, indent=4)    
\ No newline at end of file