power law exponent heat map in the uncontrolled growth phase

basic_model_only_density/rho0_tmax50/heatmap_powerlawexp.py

+#!/usr/bin/env python3
+
+import dolfin as df
+import numpy as np
+import glob
+import json
+import h5py
+import progressbar
+import os
+import argparse
+import matplotlib.pyplot as plt
+from matplotlib.backend_bases import MouseButton
+from scipy import spatial
+
+def length_heatmap(DIR):
+    J = glob.glob(DIR+'/*.json') # list of files with this path 
+
+    # array for param1, param2, quantity of interest
+    min_peclet = 0.0
+    max_peclet = 3.5
+    peclet_step = 0.5
+
+    min_rt = 0.5
+    max_rt = 5.0
+    rt_step = 0.5
+
+    peclet = np.arange(min_peclet, max_peclet + peclet_step, peclet_step)
+    retardationtime = np.arange(min_rt, max_rt + rt_step, rt_step)
+    powerlaw_exponent = np.zeros((len(peclet), len(retardationtime)))
+
+    # check data files
+    for j in J:
+        timestamp = os.path.basename(j).split('_')[0]
+        dataFile = ['%s_density.xdmf' % timestamp,
+                '%s_density.h5' % timestamp,
+                '%s_velocity.xdmf' % timestamp,
+                '%s_velocity.h5' % timestamp,
+                '%s_displacement.xdmf' % timestamp,
+                '%s_displacement.h5' % timestamp]
+        for d in dataFile:
+            assert os.path.isfile(os.path.join(DIR,d)), '%s not found' % d
+
+        # load the json file
+        with open(j) as jFile:
+            p = json.load(jFile)
+
+            # find the peclet, rt in the file
+            pe = (p['lamda'] / (p['viscosity']*p['turnover_rho']))
+
+            rt = (p['elasticity'] / 
+                    (p['viscosity']*p['turnover_rho']))
+            
+            # extract the index of this in the arrays
+            index_pe = np.where(peclet == pe)[0]
+            index_rt = np.where(retardationtime == rt)[0]
+            
+            # load the length data
+            length = np.load('%s_length.npy' % p['timestamp'])
+
+            # time array
+            savesteps = int(p['maxtime']/p['savetime'])
+            times = np.arange(savesteps + 1) * p['savetime']
+
+            # find the quantity of interest at index found
+            if(length[-1]-length[-3]<1e-7):
+                print('controlled_growth')
+                powerlaw_exponent[index_pe, index_rt] = 0.0
+            else:
+                print('uncontrolled_growth')
+
+                # fit last 80 % of data
+                times_clipped = times[-int(8*savesteps/10):]
+                length_clipped = length[-int(8*savesteps/10):]
+                logt = np.log(times_clipped)
+                loglength = np.log(length_clipped)
+
+                # fit the data
+                m, c = np.polyfit(logt, loglength, 1)
+                powerlaw_exponent[index_pe, index_rt] = m
+
+                # save a plot of the data
+                plt.figure()
+                plt.grid()
+                plt.xlabel(r"$\log t$")
+                plt.ylabel(r"$\log L(t)$")
+                plt.loglog(times_clipped, length_clipped,'o',ms=3)
+                plt.loglog(times_clipped, np.exp(m*logt + c), color='r',linewidth = 2.5)
+                plt.savefig('%s_power_law_fit_length.png' % p['timestamp'])
+                plt.close()
+               
+    plt.figure()
+    plt.xlabel('Retardation time/birth-death time')
+    plt.ylabel('Peclet')
+    fig1 = plt.imshow(powerlaw_exponent, cmap='Greens', 
+                   extent = [np.min(retardationtime)-rt_step/2, np.max(retardationtime)+rt_step/2, 
+                             np.min(peclet)-peclet_step/2, np.max(peclet)+peclet_step/2],
+                    interpolation ='nearest', origin ='lower') 
+    plt.colorbar(fig1)
+    plt.title('Power law exponent in the uncontrolled phase for $L_i = %s$' % p['system_size'])
+    plt.show()
+
+if __name__=="__main__":
+    import argparse
+    parser = argparse.ArgumentParser()
+    parser.add_argument('-d','--directory', 
+            help='directory with json files', required=True)
+    args = parser.parse_args()
+    assert os.path.isdir(args.directory), \
+            '%s directory not found' % args.directory
+    length_heatmap(args.directory)
+