fixed linear growth

growing_domain/diffusion_on_growing_domain.py

@@ -13,22 +13,27 @@ class GrowthDiffusion(object):
         for key in parameters:
         for key in parameters:
             setattr(self, key, parameters[key])
             setattr(self, key, parameters[key])
             
             
-        # choose growth model
-        if self.growth=='none':
-            self.sigma = df.Constant(0.0)
-        elif self.growth=='linear':
-            self.sigma = df.Expression('b/(L0 + b*t)', 
-                            b=self.growth_parameter, L0=self.system_size, degree=1)
-        elif self.growth=='exponential':
-            self.sigma = df.Constant(self.growth_parameter)
 
 
         # set up mesh
         # set up mesh
         # and define the growth velocity field
         # and define the growth velocity field
         if self.dimension==1:
         if self.dimension==1:
-            self.velocity = df.Expression(('s*x[0]',), s=self.sigma, degree=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)
             self.mesh = df.IntervalMesh(self.resolution, 0, self.system_size)
+
         elif self.dimension==2:
         elif self.dimension==2:
-            self.velocity = df.Expression(('s*x[0]', 's*x[1]'), s=self.sigma, degree=1)
+            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, cells = meshzoo.disk(self.resolution, 2*self.resolution)
             points *= self.system_size
             points *= self.system_size
             self.mesh = df.Mesh()
             self.mesh = df.Mesh()
@@ -49,7 +54,7 @@ class GrowthDiffusion(object):
         self.c0 = df.Function(self.SFS) 
         self.c0 = df.Function(self.SFS) 
         tc = df.TestFunction(self.SFS)
         tc = df.TestFunction(self.SFS)
         
         
-        self.velocity.t = 0
+        # self.velocity.t = 0
         self.vel = df.project(self.velocity, self.VFS)
         self.vel = df.project(self.velocity, self.VFS)
         self.form = ( df.inner((self.c-self.c0)/self.timestep, tc) 
         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.vel*self.c0), tc) 
@@ -58,7 +63,12 @@ class GrowthDiffusion(object):
 
 
     def solve(self):
     def solve(self):
         times = np.arange(0, self.maxtime+self.timestep, self.timestep)
         times = np.arange(0, self.maxtime+self.timestep, self.timestep)
-        cFile = df.XDMFFile('concentration.xdmf')
+        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')
         
         
         # initial condition
         # initial condition
         if self.dimension==1:
         if self.dimension==1:
@@ -69,7 +79,13 @@ class GrowthDiffusion(object):
                                pi=np.pi, m=2, L=self.system_size, degree=1)
                                pi=np.pi, m=2, L=self.system_size, degree=1)
         self.c0.assign(df.project(c0, self.SFS))
         self.c0.assign(df.project(c0, self.SFS))
         # save data
         # save data
-        cFile.write_checkpoint(self.c0, 'concentration', 0.0)
+        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)
+
         # time stepping
         # time stepping
         for i in progressbar.progressbar(range(1, len(times))):
         for i in progressbar.progressbar(range(1, len(times))):
             # get velocity
             # get velocity
@@ -80,7 +96,12 @@ class GrowthDiffusion(object):
             # update
             # update
             self.c0.assign(self.c)
             self.c0.assign(self.c)
             # save data
             # save data
-            cFile.write_checkpoint(self.c0, 'concentration', times[i], append=True)
+            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)
             # move mesh
             # move mesh
             displacement = df.project(self.vel*self.timestep, self.VFS)
             displacement = df.project(self.vel*self.timestep, self.VFS)
             df.ALE.move(self.mesh, displacement)
             df.ALE.move(self.mesh, displacement)