Working diffusion on growing domains in both 1D and 2D

growing_domain/diffusion_on_growing_domain.py

@@ -2,26 +2,91 @@ import dolfin as df
 import numpy as np
 import progressbar
 import os
+import meshzoo
 
 df.set_log_level(df.LogLevel.ERROR)
 df.parameters['form_compiler']['optimize'] = True
 
 class GrowthDiffusion(object):
-    def __init__(self, parameters, mDir=''):
+    def __init__(self, parameters):
         # read in parameters
         for key in parameters:
             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
+        # and define the growth velocity field
+        if self.dimension==1:
+            self.velocity = df.Expression(('s*x[0]',), s=self.sigma, degree=1)
+            self.mesh = df.IntervalMesh(self.resolution, 0, self.system_size)
+        elif self.dimension==2:
+            self.velocity = df.Expression(('s*x[0]', 's*x[1]'), s=self.sigma, degree=1)
+            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()
+
         # create mesh, function space, define function, test function
-        self.mesh = df.Mesh(os.path.join(mDir, self.meshfile))
-        self.function_space = df.FunctionSpace(self.mesh, 'P', 1)
-        self.c, self.tc = df.Function(self.function_space), df.TestFunction(self.function_space)
+        self.SFS = df.FunctionSpace(self.mesh, 'P', 1)
+        self.VFS = df.VectorFunctionSpace(self.mesh, 'P', 1)
+        self.c = df.Function(self.SFS)
+        self.c0 = df.Function(self.SFS) 
+        tc = df.TestFunction(self.SFS)
+        
+        self.velocity.t = 0
+        self.vel = df.project(self.velocity, self.VFS)
+        self.form = ( df.inner((self.c-self.c0)/self.timestep, tc) 
+                        + df.inner(df.div(self.vel*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)
+        cFile = df.XDMFFile('concentration.xdmf')
         
-    def setup_initial_conditions():
+        # 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)
+        self.c0.assign(df.project(c0, self.SFS))
+        # save data
+        cFile.write_checkpoint(self.c0, 'concentration', 0.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))
+            # solve
+            df.solve(self.form == 0, self.c)
+            # update
+            self.c0.assign(self.c)
+            # save data
+            cFile.write_checkpoint(self.c0, 'concentration', times[i], append=True)
+            # move mesh
+            displacement = df.project(self.vel*self.timestep, self.VFS)
+            df.ALE.move(self.mesh, displacement)
 
-    def setup_weak_forms():
+        cFile.close()
     
-    def solve():
 
 if __name__ == '__main__':
     import json

growing_domain/parameters.json

 {
-    "dimension" : 1,
-    "mesh" : 32,
-    "system_size_by_2PI" : 1 ,
+    "dimension" : 2,
+    "resolution" : 5,
+    "system_size" : 1,
 
-    "Dc" : 1.0,
-    "growth" : "none",
-    "growth_parameter" : 1.0,
-    "reaction_rate" : 1.0,
+    "Dc" : 0.1,
+    "growth" : "exponential",
+    "growth_parameter" : 0.1,
+    "reaction_rate" : 0.0,
 
     "timestep" : 0.01,
-    "savetime" : 0.1,
-    "maxtime" : 10.0
+    "maxtime" : 1.0
 }