error vs dt for linear growth model

euler_error/linear_growth_euler_errors.py

+# moving domain diffusion advection equation
+# advection velocity v = alpha*x. domain also moves by thsi velocity
+# boundary condition: diffusive flux at the boundaries [0,1] is zero
+# initial condition c(x,0) = 1 + 0.2*cos(pi x/L) ; satisfies the boundary condition
+# exact solution: c(x,t) = exp(-alpha * t)(1 + 0.2 cos(pi * x/L0 *exp(-alpha*t)) * exp(-pi**2*D*(1-exp(-2*alpha*t)/(2*alpha*L0**2))))
+
+import numpy as np
+import dolfin as df
+import matplotlib.pyplot as plt
+import progressbar
+from scipy.optimize import curve_fit
+
+df.set_log_level(df.LogLevel.ERROR)
+df.parameters['form_compiler']['optimize'] = True
+
+def advection_diffusion(Nx, L, Nt, tmax, v, D):
+    # mesh, function space, function, test function
+    mesh = df.IntervalMesh(Nx, 0, L)
+    function_space = df.FunctionSpace(mesh, 'P', 1)
+    c, tc = df.Function(function_space), df.TestFunction(function_space)
+
+    # initial condition
+    c0 = df.interpolate(df.Expression('1 + 0.2*cos(m * pi*x[0]/L)', pi = np.pi, L = L, m = m, degree=1), function_space)
+
+    # arrays
+    times = np.linspace(0, tmax, Nt+1)
+    dt = times[1] - times[0]
+    x_array = np.zeros((Nt+1, Nx+1))
+    x_array[0] = mesh.coordinates()[:,0]
+    c_array = np.zeros((Nt+1, Nx+1))
+    c_array[0] = c0.compute_vertex_values(mesh)
+
+    # form
+    u = df.project(v, function_space)
+
+    form = (df.inner((c - c0)/dt, tc)
+            + df.inner((u * c0).dx(0), tc)
+            + D * df.inner(tc.dx(0), c.dx(0))) * df.dx(mesh)
+
+    for i in progressbar.progressbar(range(1, len(times))):
+        v.t = times[i]
+        u.assign(df.project(v, function_space))
+        df.solve(form == 0, c)
+        c_array[i] = c.compute_vertex_values(mesh)
+        c0.assign(c)
+        df.ALE.move(mesh, df.project(u*dt, function_space))
+        x_array[i] = mesh.coordinates()[:,0]
+        
+
+    return [c_array,x_array]
+
+Nx, L, tmax, D, b, m = 100, 1, 1, 0.1, 0.01, 2
+v = df.Expression('b*x[0]/(L + b*t)', b=b, L=L, t=0, degree=1)
+Nt_array = np.array([5, 10, 15, 20, 25, 30])
+# , 400, 800, 1600, 3200, 6400, 12800])
+dt_array = tmax/(Nt_array)
+error_in_c = np.zeros(len(Nt_array))
+
+for i in range(len(Nt_array)):
+    c, x = advection_diffusion(Nx, L, Nt_array[i], tmax, v, D)
+    times = np.linspace(0, tmax, Nt_array[i]+1)
+    c_exact = np.zeros((Nt_array[i]+1, Nx+1))
+    for j in range(0, Nt_array[i]+1):
+        l = L + b * times[j]
+        c_exact[j] = (L/l) * (1 + 0.2 * np.cos(m * np.pi * x[j]/l) * np.exp(-m**2 * np.pi**2 * D * times[j]/(L * l)))
+    error_in_c[i] = np.max(np.abs(c[-1] - c_exact[-1]))
+
+def linear_func(x, a, b):
+    return a * x + b
+    
+popt, pcov = curve_fit(linear_func, np.log(dt_array), np.log(error_in_c))
+print('The slope of log(error) vs log(dt) is', popt[0])
+plt.loglog(dt_array, error_in_c, 'bo')
+plt.xlabel('log(dt)')
+plt.ylabel('log(error)')
+plt.show()

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()
+        elif self.dimension==2: 
+            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()
\ No newline at end of file
+    gd.solve()
+    
+    with open(params['timestamp'] + '_parmeters.json', "w") as fp:
+        json.dump(params, fp, indent=4)    
\ No newline at end of file