working 2D exact soolutions with visualization, 1D visualization to be done

growing_domain/diffusion_on_growing_domain_exact_solutions.py

@@ -22,9 +22,6 @@ offscreen=(not args.onscreen)
 # parameters
 # parameters
 Nt = int(params['maxtime']/params['timestep'])
 Nt = int(params['maxtime']/params['timestep'])
 times = np.linspace(0, params['maxtime'], Nt+1)
 times = np.linspace(0, params['maxtime'], Nt+1)
-
-mesh = df.Mesh('disk.xml.gz')
-Nv = mesh.num_vertices()
 d = params['dimension']
 d = params['dimension']
 L = params['system_size']
 L = params['system_size']
 alpha = params['growth_parameter']
 alpha = params['growth_parameter']
@@ -32,35 +29,44 @@ dt = params['timestep']
 Dc = params['Dc']
 Dc = params['Dc']
 k = params['reaction_rate']
 k = params['reaction_rate']
 growth = params['growth']
 growth = params['growth']
+m = 2     # mth mode of cosine in the initial condition for d=1
+p = 1     # pth zero of bessel 1 in the initial condition for d=2
 
 
-# topology and geometry arrays
-topology_array = mesh.cells()
+# mesh
+if d==2:
+    mesh = df.Mesh('disk.xml.gz')
+else:
+    mesh = df.IntervalMesh(params['resolution'], 0, L)
+
+Nv = mesh.num_vertices()            # why is the number of vertices in the 2D mesh less than the number of cells?
 
 
+# initialize topology and geometry arrays
+topology_array = mesh.cells()
 geometry_array = np.zeros(((Nt + 1, Nv, d)))
 geometry_array = np.zeros(((Nt + 1, Nv, d)))
-x = np.reshape(mesh.coordinates()[:,0], (Nv, 1))
-y = np.reshape(mesh.coordinates()[:,1], (Nv, 1))
-geometry_array[0] = np.concatenate((x, y), axis=1)
 
 
-# r and solution array
+# initialize r and sol arrays
 r_array = np.zeros((Nt + 1, Nv))
 r_array = np.zeros((Nt + 1, Nv))
-r_array[0] = np.sqrt(mesh.coordinates()[:,0]**2 + mesh.coordinates()[:,1]**2)
 sol = np.zeros((len(times), len(r_array[0])))
 sol = np.zeros((len(times), len(r_array[0])))
 
 
-# velocity
+# velocity to move the mesh
 VFS = df.VectorFunctionSpace(mesh, 'P', 1)
 VFS = df.VectorFunctionSpace(mesh, 'P', 1)
 sigma = {'none': '0.0',
 sigma = {'none': '0.0',
         'linear': 'alpha/(L0 + alpha*t)',
         'linear': 'alpha/(L0 + alpha*t)',
         'exponential': 'alpha'}
         'exponential': 'alpha'}
 
 
 sigma = df.Expression(sigma[growth], L0 = L, alpha = alpha, t = 0, degree=1)
 sigma = df.Expression(sigma[growth], L0 = L, alpha = alpha, t = 0, degree=1)
-growth_direction = ('x[0]', 'x[1]')
+growth_direction = tuple(['x['+str(i)+']' for i in range(0, d)])
 growth_direction = df.Expression(growth_direction, degree=1)
 growth_direction = df.Expression(growth_direction, degree=1)
-velocity = df.project(sigma*growth_direction, VFS)
+velocity = df.project(sigma * growth_direction, VFS)
+
+# modes enetering the solution
+if d == 1:
+    mode = m * np.pi
+else:
+    mode = sc.special.jn_zeros(1, p)
 
 
 # time loop
 # time loop
 t = 0
 t = 0
-first_zero_of_bessel_1 = sc.special.jn_zeros(1, 1)
-
 for steps in range(0, Nt+1):
 for steps in range(0, Nt+1):
     # update sigma
     # update sigma
     sigma.t = t
     sigma.t = t
@@ -68,19 +74,34 @@ for steps in range(0, Nt+1):
     # update velocity
     # update velocity
     velocity.assign(df.project(sigma*growth_direction, VFS))
     velocity.assign(df.project(sigma*growth_direction, VFS))
 
 
-    # find r, geometry on solution on the mesh
-    x = np.reshape(mesh.coordinates()[:,0], (mesh.num_vertices(), 1))
-    y = np.reshape(mesh.coordinates()[:,1], (mesh.num_vertices(), 1))
-    r_array[steps] = np.sqrt(mesh.coordinates()[:,0]**2 + mesh.coordinates()[:,1]**2)
-    
-    geometry_array[steps] = np.concatenate((x,y),axis=1)
-
-    time_part = {'none': np.exp(-first_zero_of_bessel_1**2 * Dc * t/L**2 + k*t),
-        'exponential': np.exp(-first_zero_of_bessel_1**2 * Dc * (1 - np.exp(-2*alpha*t))/(2*alpha*L**2) + (k-2*alpha)*t),
-        'linear': (L/(L + alpha*t))**2 * np.exp(-first_zero_of_bessel_1**2 * Dc * t/(L*(L + alpha * t)) + k*t)}
-        
-    sol[steps] = (sc.special.j0(first_zero_of_bessel_1 * r_array[steps]/ L) 
-            * time_part[growth])
+    # find r, geometry and sol arrays on the mesh
+
+    # geometry array
+    if d==1:
+        x = np.reshape(mesh.coordinates()[:], (Nv, 1))
+        geometry_array[steps] = x
+    else:
+        x = np.reshape(mesh.coordinates()[:,0], (Nv, 1))
+        y = np.reshape(mesh.coordinates()[:,1], (Nv, 1))
+        geometry_array[steps] = np.concatenate((x,y),axis=1)
+
+    # r array
+    r_array[0] = 0
+    for j in range(0,d):
+        r_array[steps] += mesh.coordinates()[:,j]**2
+    r_array[steps] = np.sqrt(r_array[steps])
+
+    # sol array
+    time_part = {'none': np.exp(-mode**2 * Dc * t/L**2 + k*t),
+                'exponential': np.exp(-mode**2 * Dc * (1 - np.exp(-2 * alpha * t))/(2 * alpha * L**2) + (k - d * alpha) * t),
+                'linear': (L/(L + alpha * t))**d * np.exp(-mode**2 * Dc * t/(L*(L + alpha * t)) + k*t)}
+    domain_length = {'none': L,
+              'exponential': L * np.exp(alpha * t),
+              'linear': L + alpha * t}
+    if d==1:
+        sol[steps] = np.cos(mode * r_array[steps]/domain_length[growth]) * time_part[growth]
+    else:
+        sol[steps] = sc.special.j0(mode * r_array[steps]/ domain_length[growth]) * time_part[growth]
 
 
     # move the mesh
     # move the mesh
     displacement = df.project(velocity * dt, VFS)
     displacement = df.project(velocity * dt, VFS)
@@ -88,11 +109,12 @@ for steps in range(0, Nt+1):
 
 
     # update time
     # update time
     t += dt
     t += dt
-
+    
 # visualise the solution
 # visualise the solution
 n_cmap_vals = 16
 n_cmap_vals = 16
 scalar_cmap = 'viridis'
 scalar_cmap = 'viridis'
-geometry = np.dstack((geometry_array, np.zeros(geometry_array.shape[0:2])))
+geometry = np.dstack((geometry_array, np.zeros(geometry_array.shape[0:2]))) # not workin without this, check why 
+                                                                            # this is necessary
 cmin, cmax = np.min(sol), np.max(sol)
 cmin, cmax = np.min(sol), np.max(sol)
 plotter = vd.plotter.Plotter(axes=0)
 plotter = vd.plotter.Plotter(axes=0)
 poly = vd.utils.buildPolyData(geometry[0], topology_array)
 poly = vd.utils.buildPolyData(geometry[0], topology_array)