extend code + visualising dldt

basic_model_only_density/parameters.json

 {
     "resolution": 1000,
     "system_size": 1.0,
-    "elasticity": 1,
+    "elasticity": 0.5,
     "viscosity": 1,
     "zero_rho_boundary": true,
     "friction": 1.0,
-    "lamda": 1,
+    "lamda": 1.5,
     "diffusion_rho": 0.1,
     "turnover_rho": 1,
     "active_death": false,
@@ -15,5 +15,5 @@
     "noise_level": 0.0,
     "timestep": 0.01,
     "savetime": 0.1,
-    "maxtime": 1.0
+    "maxtime": 600.0
 }
\ No newline at end of file

basic_model_only_density/run_tissue.py

@@ -38,17 +38,12 @@ else:
     mesh = df.IntervalMesh(parameters['resolution'], 
                             - parameters['system_size']/2, 
                               parameters['system_size']/2)
-    print('Before reading geometry the mesh:', mesh.coordinates()[:,0])
     
     # read geometry from h5 file
     var = 'density'
     h5 = h5py.File( '%s_%s.h5' % (parameters['timestamp'], var,), 'a')
-    # print("Name of the file: ", h5.filename)
-    # print('Details of the h5 file:', [(i,j) for  i,j in h5.attrs.items()])
-    print(f'var {var}, savesteps {savesteps}')
     geometry = np.array(h5['%s/%s_%d/mesh/geometry'%(var,var, savesteps)])
     mesh.coordinates()[:,0] = geometry[:,0]
-    print('After reading geometry the mesh:',mesh.coordinates()[:,0])
 
     # Read data
     SFS = df.FunctionSpace(mesh, 'P', 1)
@@ -62,7 +57,6 @@ else:
     rhoFile.close()
 
     initTime = parameters['maxtime']
-    print('old max time=', initTime)
     parameters['maxtime'] = args.time
 
 tissue = Tissue(parameters, mesh)
@@ -74,5 +68,5 @@ if extend:
 with open(parametersFile, "w") as jsonFile:
     json.dump(parameters, jsonFile, indent=4, sort_keys=True)
 
-# from viz_tissue import visualize
-# visualize(parameters, DIR='')
\ No newline at end of file
+from viz_tissue import visualize
+visualize(parameters, DIR='')
\ No newline at end of file

basic_model_only_density/tissue.py

@@ -17,7 +17,6 @@ class Tissue(object):
                                     self.system_size/2)
         else:
             self.mesh = mesh
-            print('Mesh passed as initial condition:', self.mesh.coordinates()[:,0])
 
         scalar_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1)
         vector_element = df.VectorElement('P', self.mesh.ufl_cell(), 1)
@@ -156,7 +155,6 @@ class Tissue(object):
             self.rhoFile.write_checkpoint(rho, 'density', self.time)
 
         for steps in progressbar.progressbar(range(1, maxsteps + 1)):
-            print(self.time)
             df.solve(self.form == 0, self.function, self.bc)
             self.function0.assign(self.function)
             self.time += self.timestep            
@@ -170,17 +168,13 @@ class Tissue(object):
                                         self.time, append=True)
                 self.rhoFile.write_checkpoint(rho, 'density', 
                                             self.time, append=True)
+                
+
             # move mesh
             dr = df.project(v*self.timestep, self.function_space.sub(0).collapse())
             df.ALE.move(self.mesh, dr)
 
-            print('After solving till new maxtime, the mesh coordinates:',self.mesh.coordinates()[:,0])
-
-
-        self.uFile = None
-        self.vFile = None
-        self.rhoFile = None
-        # self.uFile.closed
-        # self.vFile.closed
-        # self.rhoFile.closed        
+        self.uFile.close()
+        self.vFile.close()
+        self.rhoFile.close()
     

basic_model_only_density/viz_tissue.py

@@ -10,7 +10,6 @@ from tempfile import TemporaryDirectory
 def visualize(params, DIR=''):
     savesteps = int(params['maxtime']/params['savetime'])
     times = np.arange(savesteps+1) * params['savetime']
-    print(times)
     
     # Read mesh geometry from h5 file
     var = 'density'
@@ -21,7 +20,6 @@ def visualize(params, DIR=''):
     topology = np.array(h5['%s/%s_0/mesh/topology'%(var,var)])
     geometry = []
     for i in progressbar.progressbar(range(len(times))):
-        print(times[i])
         geometry.append(np.array(h5['%s/%s_%d/mesh/geometry'%(var,var,i)]))
     h5.close()
     geometry = np.array(geometry)
@@ -114,7 +112,7 @@ def visualize(params, DIR=''):
             fc='#999999')
     slider.drawon = False
     slider.on_changed(update)
-    # plt.show()  
+
     print('Saving movie-...')
     FPS = 50
     movFile = os.path.join(DIR, '%s_fields.mov' % params['timestamp'])
@@ -133,7 +131,7 @@ def visualize(params, DIR=''):
     tmp_dir.cleanup()
 
     # plotting the stresses
-    fig2, axes1 = plt.subplots(3, 1, sharex=True, figsize=(8,8))
+    figure, axes1 = plt.subplots(3, 1, sharex=True, figsize=(8,8))
     axes1[-1].set_xlabel(r'$x$')
     axes1[2].set_ylabel('Elastic stress')
     axes1[1].set_ylabel('Viscous stress')
@@ -180,7 +178,7 @@ def visualize(params, DIR=''):
     for tt in progressbar.progressbar(range(len(times))):
         slider.set_val(times[tt])
         fr = get_filename("%03d.png" % tt)
-        fig2.savefig(fr, facecolor=fig2.get_facecolor(), dpi=100)
+        figure.savefig(fr, facecolor=figure.get_facecolor(), dpi=100)
     os.system(command + " " + str(fps)
                 + " -i " + tmp_dir.name + os.sep
                 + "%03d.png " + options + " " + movFile)
@@ -188,24 +186,50 @@ def visualize(params, DIR=''):
 
     # L(t) vs t
     length = geometry[:,-1,0]-geometry[:,0,0]
+    print(length)
     fig1, ax1 = plt.subplots(1, 1, figsize=(8, 8))
     ax1.set_xlabel('$t$')
     ax1.set_ylabel('$L(t)$')
     ax1.set_xlim(np.min(times), np.max(times))
-    ax1.set_ylim(np.min(length), np.max(length)+1)
+    ax1.set_ylim(np.min(length), np.max(length)+0.01)
     ax1.plot(times, length)
     plt.savefig("%s.png" %params['timestamp'])
-    # plt.show()     
     np.save('%s_length.npy' %params['timestamp'], length)
 
     fig2, ax2 = plt.subplots(1, 1, figsize=(8, 8))
     ax2.set_xlabel('$t$')
     ax2.set_ylabel('$L(t)$')
-    # ax2.set_xlim(np.min(times), np.max(times))
-    # ax2.set_ylim(np.min(length), np.max(length)+1)
     ax2.loglog(times, length)
     plt.savefig("%s_loglog.png" %params['timestamp'])
-    # plt.show()     
+
+    fig3, ax3 = plt.subplots(1, 1, figsize=(8, 8))
+    ax3.set_xlabel('$t$')
+    ax3.set_ylabel('$L(t)$')
+    ax3.semilogy(times, length)
+    plt.savefig("%s_semilog.png" %params['timestamp'])
+
+    # dldt vs t
+    dldt = np.gradient(length, times)
+    print(dldt)
+    fig11, ax11 = plt.subplots(1, 1, figsize=(8, 8))
+    ax11.set_xlabel('$t$')
+    ax11.set_ylabel('$dL/dt$')
+    ax11.set_xlim(np.min(times), np.max(times))
+    ax11.set_ylim(np.min(dldt), np.max(dldt)+0.01)
+    ax11.plot(times, dldt)
+    plt.savefig("%s_dldt.png" %params['timestamp'])
+
+    fig21, ax21 = plt.subplots(1, 1, figsize=(8, 8))
+    ax21.set_xlabel('$t$')
+    ax21.set_ylabel('$dL/dt$')
+    ax21.loglog(times, dldt)
+    plt.savefig("%s_loglog_dldt.png" %params['timestamp'])
+
+    fig31, ax31 = plt.subplots(1, 1, figsize=(8, 8))
+    ax31.set_xlabel('$t$')
+    ax31.set_ylabel('$dL/dt$')
+    ax31.semilogy(times, dldt)
+    plt.savefig("%s_semilog_dldt.png" %params['timestamp'])
 
 if __name__ == '__main__':  
     import argparse, json