works in 1D, 2D -- visualize both analytical and numerical solutions

growing_domain/diffusion_on_growing_domain_exact_solutions.py → growing_domain/diffusion_growing_domain_an_vs_num.py

@@ -4,7 +4,6 @@ import numpy as np
 import dolfin as df
 import h5py
 from matplotlib.widgets import Slider
-
 import os
 import argparse, json   
 from tempfile import TemporaryDirectory
@@ -22,33 +21,72 @@ 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)
+tmax = params['maxtime']
 d = params['dimension']
 L = params['system_size']
 alpha = params['growth_parameter']
 dt = params['timestep']
+dts = params['savetime']
 Dc = params['Dc']
 k = params['reaction_rate']
 growth = params['growth']
+timestamp = params['timestamp']
 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
 
-# mesh
-if d==2:
-    mesh = df.Mesh('disk.xml.gz')
-else:
-    mesh = df.IntervalMesh(params['resolution'], 0, L)
-
-Nv = mesh.num_vertices()         # gives params['resolution'] + 1  
-
-# topology array, initialize geometry arrays
-topology_array = mesh.cells()
-geometry_array = np.zeros(((Nt + 1, Nv, d)))
+#################### get the numerical solution, the mesh, its geometry and topology ##################################
+def get_data(params, DIR=''):
+    savesteps = int(params['maxtime']/params['savetime'])
+    times = np.arange(savesteps+1) * params['savetime']    
+
+    # Read mesh geometry from h5 file
+    var = 'concentration'
+    h5 = h5py.File(os.path.join(DIR, '%s_%s.h5' % (timestamp, var)), "r")
+    # should be in the loop if remeshing
+    topology = np.array(h5['%s/%s_0/mesh/topology'%(var,var)])
+    geometry = []
+    for i in range(len(times)):
+        geometry.append(np.array(h5['%s/%s_%d/mesh/geometry'%(var,var,i)]))
+    h5.close()
+    geometry = np.array(geometry)
+    if d==1:
+        geometry, zeros= np.dsplit(geometry, 2)
 
-# initialize r and sol arrays
-r_array = np.zeros((Nt + 1, Nv))
-sol = np.zeros((len(times), len(r_array[0])))
+    # create a mesh
+    if params['dimension']==1:
+        mesh = df.IntervalMesh(params['resolution'], 0, L)
+    else:
+        mesh = df.Mesh()
+        editor = df.MeshEditor()
+        editor.open(mesh, 
+                    type='triangle' if d==2 else 'tetrahedron', 
+                    tdim=d, gdim=d)
+        editor.init_vertices(len(geometry[0]))
+        editor.init_cells(len(topology))
+        for i in range(len(geometry[0])):
+            editor.add_vertex(i, geometry[0][i])  # vertex 0 --> (x0, y0, z0), vertex 1 --> (x1, y1, z1), ...
+        for i in range(len(topology)):
+            editor.add_cell(i, topology[i])       # cell 0 --> vertices(0,1,2), cell 1 --> ()...
+        editor.close()
+    
+    # Read concentration data
+    concentration = np.zeros((len(times), len(geometry[0])))
+    cFile = os.path.join(DIR, '%s_%s.xdmf' % (timestamp,var))
+    with df.XDMFFile(cFile) as cFile:
+        for steps in range(savesteps+1):
+            mesh.coordinates()[:] = geometry[steps]
+            VFS = df.FunctionSpace(mesh, 'P', 1)
+            c = df.Function(VFS)
+            cFile.read_checkpoint(c, var, steps)
+            concentration[steps] = c.compute_vertex_values(mesh)
+
+    return (mesh, times, geometry, topology, concentration)
+
+(mesh, times, geometry, topology, concentration) = get_data(params, DIR)
+
+#################################### get the analytical solution #########################################################
+radius = np.zeros((len(times), mesh.num_vertices()))
+sol = np.zeros_like(radius)
 
 # velocity to move the mesh
 VFS = df.VectorFunctionSpace(mesh, 'P', 1)
@@ -61,7 +99,7 @@ growth_direction = tuple(['x['+str(i)+']' for i in range(0, d)])
 growth_direction = df.Expression(growth_direction, degree=1)
 velocity = df.project(sigma * growth_direction, VFS)
 
-# modes enetering the solution
+# modes entering the solution
 if d == 1:
     mode = m * np.pi
 else:
@@ -69,29 +107,16 @@ else:
 
 # time loop
 t = 0
-for steps in range(0, Nt+1):
+
+for steps in range(0, len(times)):
     # update sigma
     sigma.t = t
 
     # update velocity
     velocity.assign(df.project(sigma*growth_direction, VFS))
 
-    # 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[steps] = 0
-    for j in range(0,d):
-        r_array[steps] += mesh.coordinates()[:,j]**2
-    r_array[steps] = np.sqrt(r_array[steps])
+    # r_array
+    radius[steps] = np.linalg.norm(geometry[steps], axis=1)
 
     # sol array
     mode_indep_time_part = {'none': np.exp(k*t),
@@ -108,32 +133,34 @@ for steps in range(0, Nt+1):
                             'exponential': L * np.exp(alpha * t),
                             'linear': L + alpha * t}
     if d==1:
-        sol[steps] = (1 + mode_dep_time_part[growth] * np.cos(mode * r_array[steps]/domain_length[growth])) * mode_indep_time_part[growth]
+        sol[steps] = (1 + mode_dep_time_part[growth] * np.cos(mode * radius[steps]/domain_length[growth])) * mode_indep_time_part[growth]
     else:
-        sol[steps] = (1 + mode_dep_time_part[growth] * sc.special.j0(mode * r_array[steps]/ domain_length[growth])) * mode_indep_time_part[growth]
+        sol[steps] = (1 + mode_dep_time_part[growth] * sc.special.j0(mode * radius[steps]/ domain_length[growth])) * mode_indep_time_part[growth]
 
     # move the mesh
     displacement = df.project(velocity * dt, VFS)
     df.ALE.move(mesh, displacement)
 
     # update time
-    t += dt
-
+    t += dts
 
-# visualise the solution
+###################### visualise the analytical solution #######################################################
 
 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_xlim(np.min(radius), np.max(radius))
+    axc.set_ylim(min(np.min(sol), np.min(concentration)), max(np.max(sol), np.max(concentration)))
     axc.set_ylabel(r'$c(x,t)$')
-    cplot, = axc.plot(r_array[0], sol[0])
+    c_exactplot, = axc.plot(radius[0], sol[0])
+    cplot, = axc.plot(radius[0], concentration[0])
 
     def update(value):
         ti = np.abs(times-value).argmin()
-        cplot.set_ydata(sol[ti])
-        cplot.set_xdata(r_array[ti])
+        c_exactplot.set_ydata(sol[ti])
+        c_exactplot.set_xdata(radius[ti])
+        cplot.set_ydata(concentration[ti])
+        cplot.set_xdata(radius[ti])
         plt.draw()
         
     sax = plt.axes([0.1, 0.92, 0.7, 0.02])
@@ -145,17 +172,19 @@ if d==1:
     plt.show()
 
 else:
+    # for heat map
     n_cmap_vals = 16
     scalar_cmap = 'viridis'
-    geometry = np.dstack((geometry_array, np.zeros(geometry_array.shape[0:2]))) 
+    geometry = np.dstack((geometry, np.zeros(geometry.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)
+    poly = vd.utils.buildPolyData(geometry[0], topology)
     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$', 
+    scalar_actor.add_scalarbar(title = r'$c_{analytical}$', 
             pos=(0.8, 0.04), nlabels=2,
             # titleYOffset=15, titleFontSize=28, size=(100, 600)
             )
@@ -176,24 +205,48 @@ else:
 
     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()
+    fig1, axc1 = plt.subplots(1,1, figsize=(8,8))
+    axc1.set_xlabel(r'$r$')
+    axc1.set_xlim(np.min(radius), np.max(radius))
+    axc1.set_ylim(np.min(sol), np.max(sol))
+    axc1.set_ylabel(r'$c(r,t)$')
+    c_exactplot, = axc1.plot(radius[0], sol[0], 'o', ms=3, label='analytical solution')
+    cplot, = axc1.plot(radius[0], concentration[0], 'o', ms=3, label='numerical solution')
+    def update(value):
+        ti = np.abs(times-value).argmin()
+        c_exactplot.set_ydata(sol[ti])
+        c_exactplot.set_xdata(radius[ti])
+        cplot.set_ydata(concentration[ti])
+        cplot.set_xdata(radius[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)
+    axc1.legend()
+    plt.show()
+
+# make movie
+if offscreen:
+    FPS = 10
+    movFile = '%s_analytical.mov' % timestamp
+    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()