plots heatmap for final length, added length, relative change in length in the controlled phase

basic_model_only_density/heatmap_length.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)
+    final_length = np.zeros((len(peclet), len(retardationtime)))
+    added_length = np.zeros((len(peclet), len(retardationtime)))
+    relative_change_in_length = 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'])
+
+            # find the quantity of interest at index found
+            if(length[-1]-length[-3]<1e-7):
+                print('controlled_growth')
+                final_length[index_pe, index_rt] = length[-1]
+                added_length[index_pe, index_rt] = length[-1] - length[0]
+                relative_change_in_length[index_pe, index_rt] = (length[-1] - length[0])/length[0]
+            else:
+                print('uncontrolled_growth')
+                final_length[index_pe, index_rt] = 0.0
+                added_length[index_pe, index_rt] = 0.0
+                relative_change_in_length[index_pe, index_rt] = 0.0
+               
+
+    plt.figure()
+    plt.xlabel('Retardation time/birth-death time')
+    plt.ylabel('Peclet')
+    fig1 = plt.imshow(final_length, 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.title('$L_f$ in the controlled growth phase for $L_i = 1$')
+    plt.show()
+
+    plt.figure()
+    plt.xlabel('Retardation time/birth-death time')
+    plt.ylabel('Peclet')
+    fig2 = plt.imshow(added_length, 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.title('$\Delta L$ in the controlled growth phase for $L_i = 1$')
+    plt.show()
+
+    plt.figure()
+    plt.xlabel('Retardation time/birth-death time')
+    plt.ylabel('Peclet')
+    fig3 = plt.imshow(relative_change_in_length, 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.title('$\Delta L/L_i$ in the controlled growth phase for $L_i = 1$')
+    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)
+