working code for 1D, 2D exact solutions for diffusion on a growing domain for all growth models

growing_domain/diffusion_on_growing_domain_exact_solutions.py

@@ -3,6 +3,8 @@ import scipy as sc
 import numpy as np
 import dolfin as df
 import h5py
+from matplotlib.widgets import Slider
+
 import os
 import argparse, json   
 from tempfile import TemporaryDirectory
@@ -38,9 +40,9 @@ if d==2:
 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?
+Nv = mesh.num_vertices()         # gives params['resolution'] + 1  
 
-# initialize topology and geometry arrays
+# topology array, initialize geometry arrays
 topology_array = mesh.cells()
 geometry_array = np.zeros(((Nt + 1, Nv, d)))
 
@@ -63,7 +65,7 @@ velocity = df.project(sigma * growth_direction, VFS)
 if d == 1:
     mode = m * np.pi
 else:
-    mode = sc.special.jn_zeros(1, p)
+    mode = sc.special.jn_zeros(1, p)  # pth zero of bessel 1
 
 # time loop
 t = 0
@@ -86,22 +88,29 @@ for steps in range(0, Nt+1):
         geometry_array[steps] = np.concatenate((x,y),axis=1)
 
     # r array
-    r_array[0] = 0
+    r_array[steps] = 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)}
+    mode_indep_time_part = {'none': np.exp(k*t),
+                            'exponential': np.exp((k - d * alpha) * t),
+                            'linear': (L/(L + alpha * t))**d * np.exp(k*t)
+                           }
+
+    mode_dep_time_part   = {'none': np.exp(-mode**2 * Dc * t/L**2),
+                            'exponential': np.exp(-mode**2 * Dc * (1 - np.exp(-2 * alpha * t))/(2 * alpha * L**2)),
+                            'linear': np.exp(-mode**2 * Dc * t/(L*(L + alpha * 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]
+        sol[steps] = (1 + mode_dep_time_part[growth] * np.cos(mode * r_array[steps]/domain_length[growth])) * mode_indep_time_part[growth]
     else:
-        sol[steps] = sc.special.j0(mode * r_array[steps]/ domain_length[growth]) * time_part[growth]
+        sol[steps] = (1 + mode_dep_time_part[growth] * sc.special.j0(mode * r_array[steps]/ domain_length[growth])) * mode_indep_time_part[growth]
 
     # move the mesh
     displacement = df.project(velocity * dt, VFS)
@@ -110,40 +119,65 @@ for steps in range(0, Nt+1):
     # update time
     t += dt
 
+
 # visualise the solution
-n_cmap_vals = 16
-scalar_cmap = 'viridis'
-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)
-plotter = vd.plotter.Plotter(axes=0)
-poly = vd.utils.buildPolyData(geometry[0], topology_array)
-scalar_actor = vd.mesh.Mesh(poly)
-scalar_actor.pointdata['concentration'] = sol[0]
-scalar_actor.cmap(scalar_cmap, sol[0], vmin=cmin, vmax=cmax, n=n_cmap_vals)
-scalar_actor.add_scalarbar(title = r'$c$', 
+
+if d==1:
+    fig, axc = plt.subplots(1,1, figsize=(8,8))
+    axc.set_xlabel(r'$x$')
+    axc.set_xlim(np.min(r_array), np.max(r_array))
+    axc.set_ylim(np.min(sol), np.max(sol))
+    axc.set_ylabel(r'$c(x,t)$')
+    cplot, = axc.plot(r_array[0], sol[0])
+
+    def update(value):
+        ti = np.abs(times-value).argmin()
+        cplot.set_ydata(sol[ti])
+        cplot.set_xdata(r_array[ti])
+        plt.draw()
+        
+    sax = plt.axes([0.1, 0.92, 0.7, 0.02])
+    slider = Slider(sax, r'$t/\tau$', min(times), max(times),
+            valinit=min(times), valfmt='%3.1f',
+            fc='#999999')
+    slider.drawon = False
+    slider.on_changed(update)
+    plt.show()
+
+else:
+    n_cmap_vals = 16
+    scalar_cmap = 'viridis'
+    geometry = np.dstack((geometry_array, np.zeros(geometry_array.shape[0:2]))) 
+                                                                            
+    cmin, cmax = np.min(sol), np.max(sol)
+    plotter = vd.plotter.Plotter(axes=0)
+    poly = vd.utils.buildPolyData(geometry[0], topology_array)
+    scalar_actor = vd.mesh.Mesh(poly)
+    scalar_actor.pointdata['concentration'] = sol[0]
+    scalar_actor.cmap(scalar_cmap, sol[0], vmin=cmin, vmax=cmax, n=n_cmap_vals)
+    scalar_actor.add_scalarbar(title = r'$c$', 
             pos=(0.8, 0.04), nlabels=2,
             # titleYOffset=15, titleFontSize=28, size=(100, 600)
             )
-plotter += scalar_actor
+    plotter += scalar_actor
 
-def update(idx):
+    def update(idx):
         scalar_actor.points(pts=geometry[idx], transformed=False)
         scalar_actor.pointdata['concentration'] = sol[idx]
 
-def slider_update(widget, event):
+    def slider_update(widget, event):
         value = widget.GetRepresentation().GetValue()
         idx = (abs(times-value)).argmin()
         update(idx)
             
-slider = plotter.add_slider(slider_update, pos=[(0.1,0.94), (0.5,0.94)],
+    slider = plotter.add_slider(slider_update, pos=[(0.1,0.94), (0.5,0.94)],
                         xmin=times[0], xmax=times.max(), 
                         value=times[0], title=r"$t/\tau$")
 
-vd.show(interactive=(not offscreen), zoom=0.8)
+    vd.show(interactive=(not offscreen), zoom=0.8)
 
-# make movie
-if offscreen:
+    # make movie
+    if offscreen:
         FPS = 10
         movFile = '%s.mov' % dt
         fps = float(FPS)