working code for exact solutions in 2D for all growth models

growing_domain/diffusion_on_growing_domain_exact_solutions.py

+import matplotlib.pyplot as plt
+import scipy as sc
+import numpy as np
+import dolfin as df
+import h5py
+import os
+import argparse, json   
+from tempfile import TemporaryDirectory
+import vedo as vd
+import argparse, json
+
+parser = argparse.ArgumentParser()
+parser.add_argument('-j','--jsonfile', help='parameters file', required=True)
+parser.add_argument('--onscreen', action=argparse.BooleanOptionalAction)
+args = parser.parse_args()
+
+with open(args.jsonfile) as jsonFile:
+    params = json.load(jsonFile)
+DIR = os.path.dirname(args.jsonfile)
+offscreen=(not args.onscreen)
+
+# parameters
+Nt = int(params['maxtime']/params['timestep'])
+times = np.linspace(0, params['maxtime'], Nt+1)
+
+mesh = df.Mesh('disk.xml.gz')
+Nv = mesh.num_vertices()
+d = params['dimension']
+L = params['system_size']
+alpha = params['growth_parameter']
+dt = params['timestep']
+Dc = params['Dc']
+k = params['reaction_rate']
+growth = params['growth']
+
+# topology and geometry arrays
+topology_array = mesh.cells()
+
+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
+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])))
+
+# velocity
+VFS = df.VectorFunctionSpace(mesh, 'P', 1)
+sigma = {'none': '0.0',
+        'linear': 'alpha/(L0 + alpha*t)',
+        'exponential': 'alpha'}
+
+sigma = df.Expression(sigma[growth], L0 = L, alpha = alpha, t = 0, degree=1)
+growth_direction = ('x[0]', 'x[1]')
+growth_direction = df.Expression(growth_direction, degree=1)
+velocity = df.project(sigma*growth_direction, VFS)
+
+# time loop
+t = 0
+first_zero_of_bessel_1 = sc.special.jn_zeros(1, 1)
+
+for steps in range(0, Nt+1):
+    # update sigma
+    sigma.t = t
+
+    # update velocity
+    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])
+
+    # move the mesh
+    displacement = df.project(velocity * dt, VFS)
+    df.ALE.move(mesh, displacement)
+
+    # 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])))
+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
+
+def update(idx):
+    scalar_actor.points(pts=geometry[idx], transformed=False)
+    scalar_actor.pointdata['concentration'] = sol[idx]
+    
+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)],
+                        xmin=times[0], xmax=times.max(), 
+                        value=times[0], title=r"$t/\tau$")
+
+vd.show(interactive=(not offscreen), zoom=0.8)
+
+# make movie
+if offscreen:
+    FPS = 10
+    movFile = '%s.mov' % dt
+    fps = float(FPS)
+    command = "ffmpeg -y -r"
+    options = "-b:v 3600k -qscale:v 4 -vcodec mpeg4"
+    tmp_dir = TemporaryDirectory()
+    get_filename = lambda x: os.path.join(tmp_dir.name, x)
+    for tt in range(len(times)):
+        idx = (abs(times-times[tt])).argmin()
+        update(idx)
+        slider.GetRepresentation().SetValue(times[tt])
+        fr = get_filename("%03d.png" % tt)
+        vd.io.screenshot(fr)
+    os.system(command + " " + str(fps)
+                + " -i " + tmp_dir.name + os.sep
+                + "%03d.png " + options + " " + movFile)
+    tmp_dir.cleanup()
+    
+