Working 1D/2D/3D all growth models. Need to benchmark with exact solution

growing_domain/diffusion_on_growing_domain.py

 import dolfin as df
+import mshr as ms
 import numpy as np
 import progressbar
 import os
-import meshzoo
 
 df.set_log_level(df.LogLevel.ERROR)
 df.parameters['form_compiler']['optimize'] = True
@@ -12,40 +12,27 @@ class GrowthDiffusion(object):
         # 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)])
+        #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
-        # and define the growth velocity field
         if self.dimension==1:
-            if self.growth=='none':
-                self.velocity = df.Expression(('s*x[0]',), s=df.Constant(0), t=0, degree=0)
-            if self.growth == 'linear':
-                self.velocity = df.Expression(('b*x[0]/(L0 + b*t)',), b = self.growth_parameter, L0=self.system_size, t=0, degree=0)
-            if self.growth=='exponential':
-                self.velocity = df.Expression(('s*x[0]'), s=df.Constant(self.growth_parameter), t=0, degree=0)
-
             self.mesh = df.IntervalMesh(self.resolution, 0, self.system_size)
-
         elif self.dimension==2: 
-            if self.growth=='none':
-                self.velocity = df.Expression(('s*x[0]','s*x[1]'), s=df.Constant(0), t=0, degree=0)
-            if self.growth == 'linear':
-                self.velocity = df.Expression(('b*x[0]/(L0 + b*t)','b*x[1]/(L0 + b*t)'), b = self.growth_parameter, L0=self.system_size, t=0, degree=0)
-            if self.growth=='exponential':
-                self.velocity = df.Expression(('s*x[0]', 's*x[1]'), s=df.Constant(self.growth_parameter), t=0, degree=0)
-
-            points, cells = meshzoo.disk(self.resolution, 2*self.resolution)
-            points *= self.system_size
-            self.mesh = df.Mesh()
-            e = df.MeshEditor()
-            e.open(self.mesh, type='triangle', tdim=2, gdim=2)
-            e.init_vertices(len(points))
-            e.init_cells(len(cells))
-            for n in range(len(points)):
-                e.add_vertex(n, [points[n, 0], points[n, 1]])
-            for n in range(len(cells)):
-                e.add_cell(n, cells[n])
-            e.close()
+            geometry = ms.Circle(df.Point(0, 0), self.system_size)
+            self.mesh = ms.generate_mesh(geometry, self.resolution)     
+        elif self.dimension==3:
+            geometry = ms.Sphere(df.Point(0, 0, 0), self.system_size)
+            self.mesh = ms.generate_mesh(geometry, self.resolution)
 
         # create mesh, function space, define function, test function
         self.SFS = df.FunctionSpace(self.mesh, 'P', 1)
@@ -54,69 +41,59 @@ class GrowthDiffusion(object):
         self.c0 = df.Function(self.SFS) 
         tc = df.TestFunction(self.SFS)
         
-        # self.velocity.t = 0
-        self.vel = df.project(self.velocity, self.VFS)
+        self.velocity = df.project(self.sigma * self.growth_direction, self.VFS)
         self.form = ( df.inner((self.c-self.c0)/self.timestep, tc) 
-                        + df.inner(df.div(self.vel*self.c0), tc) 
+                        + df.inner(df.div(self.velocity*self.c0), tc) 
                         + self.Dc * df.inner(df.nabla_grad(self.c), df.nabla_grad(tc))
                         - df.inner(self.reaction_rate * self.c, tc) ) * df.dx
 
     def solve(self):
         times = np.arange(0, self.maxtime+self.timestep, self.timestep)
-        if self.growth =='linear':
-            cFile = df.XDMFFile('concentration_linear.xdmf')
-        elif self.growth == 'exponential':
-            cFile = df.XDMFFile('concentration_exponential.xdmf')
-        elif self.growth == 'none':
-            cFile = df.XDMFFile('concentration_no_growth.xdmf')
+        fname = params['timestamp'] + '_concentration' 
+        cFile = df.XDMFFile(fname+ '.xdmf')
+
         
         # initial condition
-        if self.dimension==1:
-            c0 = df.Expression('1 + 0.1*cos(pi*m*x[0]/L)', 
-                               pi=np.pi, m=1, L=self.system_size, degree=1)
-        elif self.dimension==2:
-            c0 = df.Expression('1 + 0.1*cos(pi*m*sqrt(x[0]*x[0]+x[1]*x[1])/L)', 
-                               pi=np.pi, m=2, L=self.system_size, degree=1)
+        r2 = "+".join(['x['+str(i)+']*x['+str(i)+']' for i in range(self.dimension)])
+        c0 = '1 + 0.1*cos(pi*m*(%s)/L)' % r2           
+        c0 = df.Expression(c0, pi=np.pi, m=1, L=self.system_size, degree=1)            
         self.c0.assign(df.project(c0, self.SFS))
         # save data
-        if self.growth=='linear':
-            cFile.write_checkpoint(self.c0, 'concentration_linear', 0.0)
-        elif self.growth=='exponential':
-            cFile.write_checkpoint(self.c0, 'concentration_exponential', 0.0)
-        elif self.growth == 'none':
-            cFile.write_checkpoint(self.c0, 'concentration_no_growth', 0.0)
+        cFile.write_checkpoint(self.c0, fname, times[0])
 
         # time stepping
         for i in progressbar.progressbar(range(1, len(times))):
             # get velocity
-            self.velocity.t = times[i-1]
-            self.vel.assign(df.project(self.velocity, self.VFS))
+            self.sigma.t = times[i-1]
+            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 self.growth=='linear':
-                cFile.write_checkpoint(self.c0, 'concentration_linear', times[i], append=True)
-            elif self.growth =='exponential':
-                cFile.write_checkpoint(self.c0, 'concentration_exponential', times[i], append=True)
-            elif self.growth == 'none':
-                cFile.write_checkpoint(self.c0, 'concentration_no_growth', times[i], append=True)
+            cFile.write_checkpoint(self.c0, fname, times[i], append=True)
             # move mesh
-            displacement = df.project(self.vel*self.timestep, self.VFS)
+            displacement = df.project(self.velocity*self.timestep, self.VFS)
             df.ALE.move(self.mesh, displacement)
 
         cFile.close()
     
 
 if __name__ == '__main__':
-    import json
+    import json, datetime
     assert os.path.isfile('parameters.json'), 'parameters.json file not found'
     with open('parameters.json') as jsonFile:
         params = json.load(jsonFile)
         
     # parse parameters
-    assert params['growth'] in ('none', 'linear', 'exponential'), 'Unknown growth model'
+    assert params['dimension'] in (1,2,3)
+    assert params['growth'] in ('none','linear','exponential'), 'Unknown growth model'
+    
+    timestamp = datetime.datetime.now().strftime("%d%m%y-%H%M%S")
+    params['timestamp'] = timestamp
         
     gd = GrowthDiffusion(params)
     gd.solve()
+    
+    with open(params['timestamp'] + '_parmeters.json', "w") as fp:
+        json.dump(params, fp, indent=4)    
\ No newline at end of file