Merge branch 'master' of http://gitlab.icts.res.in/jigyasa.watwani/growth-pattern-control

growing_domain/growing_domain_patterns/turing_patterns_on_a_growing_domain

-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 GrowthDiffusion(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.Mesh('disk.xml.gz')
-        elif self.dimension==3:
-            self.mesh = df.Mesh('sphere.xml.gz')
-        
-        # create mesh, function space, define function, test function
-        self.SFS = df.FunctionSpace(self.mesh, 'P', 1)
-        self.VFS = df.VectorFunctionSpace(self.mesh, 'P', 1)
-        self.c1 = df.Function(self.SFS)
-        self.c2 = df.Function(self.SFS)
-        self.c10 = df.Function(self.SFS) 
-        self.c20 = df.Function(self.SFS) 
-        tc1 = df.TestFunction(self.SFS)
-        tc2 = df.TestFunction(self.SFS)
-
-        def diffusion(conc, tconc):
-            return df.inner(df.nabla_grad(conc), df.nabla_grad(tconc))
-
-        def advection(velocity, conc0, tconc):
-            return df.inner(df.div(velocity * conc0), tconc)
-        
-        def reaction_c1(a, b, c1, c2, tconc):
-            return df.inner(a - (b+1)*c1 + c1**2*c2, tconc)
-
-        def reaction_c2(a, b, c1, c2, tconc):
-            return df.inner(b*c1 - c1**2*c2, tconc)
-        
-        self.velocity = df.project(self.sigma * self.growth_direction, self.VFS)
-
-        form1 = df.inner((self.c1 - self.c10)/self.timestep, tc1) 
-                        + advection(self.velocity, self.c10, tc1)
-                        + self.Dc1 * diffusion(self.c1, tc1)
-                        - reaction_c1(self.a, self.b, self.c1, self.c2, tc1))
-
-   
-        form2 = df.inner((self.c2 - self.c20)/self.timestep, tc2) 
-                        + advection(self.velocity, self.c20, tc2)
-                        + self.Dc2 * diffusion(self.c2, tc2)
-                        - reaction_c2(self.a, self.b, self.c1, self.c2, tc1))   
-
-        self.form = ( form1 + form2 ) * df.dx
-
-
-    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' %params['timestamp'])
-
-        # initial condition
-        r2 = "+".join(['x['+str(i)+']*x['+str(i)+']' for i in range(self.dimension)])
-        if self.dimension==2:
-            c10 = '1 + jn(0, m*sqrt(%s)/L)' % r2           
-            c10 = df.Expression(c10, m=sc.special.jn_zeros(1,1), L=self.system_size, degree=1)
-            c20 = '1 + jn(0, m*sqrt(%s)/L)' % r2           
-            c20 = df.Expression(c20, m=sc.special.jn_zeros(1,2), L=self.system_size, degree=1)
-        else:
-            c10 = '1 + cos(m*pi*sqrt(%s)/L)' % r2  
-            c10 = df.Expression(c10, pi=np.pi, m=2, L=self.system_size, degree=1)
-            c20 = '1 + cos(m*pi*sqrt(%s)/L)' % r2  
-            c20 = df.Expression(c20, pi=np.pi, m=3, L=self.system_size, degree=1)
-        self.c10.assign(df.project(c10, self.SFS))
-        self.c20.assign(df.project(c20, self.SFS))
-
-        # save data
-        c1File.write_checkpoint(self.c10, 'c1', self.time)
-        c2File.write_checkpoint(self.c20, 'c2', self.time)
-
-        # time stepping
-        for steps in 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.c)
-            # update
-            self.c0.assign(self.c)
-            # save data
-            if steps % savesteps == 0:
-                cFile.write_checkpoint(self.c0, 'concentration', 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
-        cFile.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,3)
-    assert params['growth'] in ('none','linear','exponential'), 'Unknown growth model'
-    if params['dimension'] in (2,3):
-        del params['resolution']  
-    
-    timestamp = datetime.datetime.now().strftime("%d%m%y-%H%M%S")
-    params['timestamp'] = timestamp
-        
-    gd = GrowthDiffusion(params)
-    gd.solve()
-    
-    with open(params['timestamp'] + '_parameters.json', "w") as fp:
-        json.dump(params, fp, indent=4)    
\ No newline at end of file