rewritten with only stress and density, substituted velocity in terms of stress

basic_model_only_density/tissue_maxwell.py

-import numpy as np
+zimport numpy as np
 import dolfin as df
 import progressbar
 import os
@@ -25,114 +25,69 @@ class Tissue(object):
 
         scalar_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1)
 
-        # v, rho, stress
-        mixed_element = df.MixedElement([scalar_element, 
-                                         scalar_element,
-                                         scalar_element])
+        # stress, rho
+        mixed_element = df.MixedElement([scalar_element, scalar_element])
 
         # define function space with this mixed element
-        self.function_space = df.FunctionSpace(self.mesh, mixed_element)
         
-        if self.zero_rho_boundary:
-            self.rho_boundary = 0.0
-        else:
-            self.rho_boundary = self.average_rho
+        self.function_space = df.FunctionSpace(self.mesh, mixed_element)
 
-        self.bc = df.DirichletBC(self.function_space.sub(1), 
-                                 df.Constant(self.rho_boundary), 
-                                 'on_boundary')
+        bc1 = df.DirichletBC(self.function_space.sub(0), df.Constant(0.0), 'on_boundary')
+        bc2 = df.DirichletBC(self.function_space.sub(1), df.Constant(0.0), 'on_boundary')
+        self.bc = ([bc1, bc2])
           
         self.function0 = df.Function(self.function_space)
         self.function  = df.Function(self.function_space)
 
-    def active_stress(self, rho):
-        if self.active_death:
-            return (-self.lamda * rho * (rho - self.active_stress_setpoint)/(rho + self.saturation_rho) 
-                   )
-        else:
-            return (-self.lamda * rho /(rho + self.saturation_rho) 
-                   )
-   
-    def fluid_stress(self, v):
-        return self.viscosity * v.dx(0)
-    
-    def diffusion_reaction_rho(self, rho, trho):
-        return (self.diffusion_rho * df.inner(rho.dx(0), trho.dx(0)) 
-               - self.turnover_rho * df.inner(rho * (1 - rho/self.average_rho),
-                                               trho)
-                )
-
     def setup_initial_conditions(self):
       
-        # ic for rho
-        if self.zero_rho_boundary:
-            base_rho = 0.0
-        else:
-            base_rho = self.average_rho
+       # ic for stress
+        stress0 = df.interpolate(df.Expression(
+                            'cos(PI*x[0]/L)*cos(PI*x[0]/L)', 
+                            L=self.system_size,
+                            PI=np.pi, degree=1), 
+                            self.function_space.sub(0).collapse())
         
         rho0 = df.interpolate(df.Expression(
-                            'base_rho + cos(PI*x[0]/L)*cos(PI*x[0]/L)', 
+                            'cos(PI*x[0]/L)*cos(PI*x[0]/L)', 
                             L=self.system_size,
-                            base_rho=base_rho,
-                            degree=1, PI=np.pi), 
-                            self.function_space.sub(1).collapse())
-
-        # ic for stress
-        stress0 = df.Constant(0.0)
-        stress0 = df.interpolate(stress0, self.function_space.sub(2).collapse())
-
-        # add noise
-        noise_rho = (self.noise_level 
-                     * (2 * np.random.random(rho0.vector().size()) - 1))
-        noise_stress = (self.noise_level 
-                * (2 * np.random.random(stress0.vector().size()) - 1))
-        
-        rho0.vector()[:] = rho0.vector()[:] + noise_rho
-        stress0.vector()[:] = stress0.vector()[:] + noise_stress
-
-        # find velocity        
-        VFS = self.function_space.sub(0).collapse()
-        v0  = df.Function(VFS)
-        tv0 = df.TestFunction(VFS)
-        v0form = (self.friction * df.inner(v0, tv0) 
-                 + df.inner(stress0, tv0.dx(0))
-                 ) * df.dx
+                            PI=np.pi, degree=1), 
+                            self.function_space.sub(1).collapse()
+                            )
+        df.assign(self.function0, [stress0, rho0])
 
-        df.solve(v0form == 0, v0)
-        df.assign(self.function0, [v0, rho0, stress0])
+    def active_stress(self, rho):
+        return -self.lamda * rho/(rho + self.saturation_rho)
 
     def setup_weak_forms(self):
-        v0, rho0, stress0  = df.split(self.function0)
-        v, rho, stress  = df.split(self.function)
-        tv, trho, tstress  = df.TestFunctions(self.function_space)
-                  
-        vform = (self.friction * df.inner(v, tv) 
-                + df.inner(stress, tv.dx(0))
+        stress0, rho0 = df.split(self.function0)
+        stress, rho = df.split(self.function)
+        
+        tstress, trho  = df.TestFunctions(self.function_space)
+                  
+        stressform = ( df.inner(stress - self.active_stress(rho), tstress)
+                      + (self.viscosity/self.elasticity) *
+                      df.inner((1/self.timestep) * (stress - stress0 - self.active_stress(rho) + self.active_stress(rho0)), tstress)
+                     + (self.viscosity/(self.elasticity*self.friction)) *
+                      df.inner(stress0.dx(0) * (stress0.dx(0) - self.active_stress(rho0).dx(0)), tstress)
+                     + (self.viscosity/self.friction) *
+                      df.inner(stress.dx(0), tstress.dx(0))
                      ) 
         
         rhoform = (df.inner((rho - rho0)/self.timestep, trho)
-                  + df.inner((v0 * rho0).dx(0), trho)
-                  + self.diffusion_reaction_rho(rho, trho)
-                  )
-
-        stressform = ( (self.viscosity/self.elasticity) * 
-                       ( df.inner((stress - self.active_stress(rho)- stress0 + self.active_stress(rho0))/self.timestep, tstress)
-                        + df.inner(v0 * (stress0 - self.active_stress(rho0)).dx(0), tstress)
-                       )
-                        + df.inner(stress - self.active_stress(rho), tstress) 
-                        - df.inner(self.fluid_stress(v), tstress)
+                   + (1/self.friction) * df.inner((rho0 * stress0.dx(0)).dx(0), trho)
+                   + self.diffusion_rho * df.inner(rho.dx(0), trho.dx(0))
+                   - self.turnover_rho * df.inner(rho * (1 - rho/self.average_rho), trho)
                     )
         
-        self.form = (vform + rhoform + stressform) * df.dx
+        self.form = (stressform + rhoform)*df.dx
            
 
     def solve(self, DIR=''):
-        self.vFile   = df.XDMFFile(os.path.join(DIR, 
-                            '%s_velocity.xdmf' % self.timestamp))
-        self.rhoFile = df.XDMFFile(os.path.join(DIR, 
-                            '%s_density.xdmf' % self.timestamp))
         self.stressFile = df.XDMFFile(os.path.join(DIR, 
                             '%s_stress.xdmf' % self.timestamp))
+        self.rhoFile = df.XDMFFile(os.path.join(DIR, 
+                            '%s_density.xdmf' % self.timestamp))
         
         self.setup_initial_conditions()
         
@@ -142,31 +97,26 @@ class Tissue(object):
         self.time = 0.0
         savesteps = int(self.savetime/self.timestep)
         maxsteps = int(self.maxtime/self.timestep)
-        v, rho, stress = self.function0.split(deepcopy=True)
-
-        self.vFile.write_checkpoint(v, 'velocity', self.time)
-        self.rhoFile.write_checkpoint(rho, 'density', self.time)
+        stress, rho = self.function0.split(deepcopy=True)
         self.stressFile.write_checkpoint(stress, 'stress', self.time)
+        self.rhoFile.write_checkpoint(rho, 'density', self.time)
 
         for steps in progressbar.progressbar(range(1, maxsteps + 1)):
             df.solve(self.form == 0, self.function, self.bc)
             self.function0.assign(self.function)
             self.time += self.timestep            
-            v, rho, stress = self.function0.split(deepcopy=True)
+            stress, rho = self.function0.split(deepcopy=True)
             if steps % savesteps == 0:
-                self.vFile.write_checkpoint(v, 'velocity', 
+                self.stressFile.write_checkpoint(stress, 'stress', 
                                             self.time, append=True)
                 self.rhoFile.write_checkpoint(rho, 'density', 
                                               self.time, append=True)
-                self.stressFile.write_checkpoint(stress, 'stress', 
-                                              self.time, append=True)
             # move mesh
-            dr = df.project(v * self.timestep, self.function_space.sub(0).collapse())
+            dr = df.project((1/self.friction)*(stress.dx(0)) * self.timestep, self.function_space.sub(0).collapse())
             df.ALE.move(self.mesh, dr)
 
-        self.vFile.close()
-        self.rhoFile.close()  
         self.stressFile.close()
+        self.rhoFile.close()        
     
 if __name__ == '__main__':
     import json, datetime

basic_model_only_density/viz_tissue_maxwell.py

@@ -30,12 +30,9 @@ def visualize(params, DIR=''):
                               params['system_size']/2)
     
     # Read data
-    v = np.zeros((len(times), mesh.num_vertices()))
     rho = np.zeros((len(times), mesh.num_vertices()))
-    stress = np.zeros_like(v)
+    stress = np.zeros_like(rho)
 
-    vFile   = df.XDMFFile(os.path.join(DIR, 
-                                    '%s_velocity.xdmf' % params['timestamp']))
     rhoFile = df.XDMFFile(os.path.join(DIR, 
                                     '%s_density.xdmf' % params['timestamp']))
     stressFile = df.XDMFFile(os.path.join(DIR, 
@@ -47,17 +44,14 @@ def visualize(params, DIR=''):
     for steps in progressbar.progressbar(range(savesteps+1)):
         mesh.coordinates()[:] = geometry[steps]
         SFS = df.FunctionSpace(mesh, 'P', 1)
-        vi, rhoi, stressi = df.Function(SFS), df.Function(SFS), df.Function(SFS)
+        rhoi, stressi = df.Function(SFS), df.Function(SFS)
             
-        vFile.read_checkpoint(vi, 'velocity', steps)
         rhoFile.read_checkpoint(rhoi, 'density', steps)
         stressFile.read_checkpoint(stressi, 'stress', steps)
 
-        v[steps] = vi.compute_vertex_values(mesh)
         rho[steps] = rhoi.compute_vertex_values(mesh)
         stress[steps] = stressi.compute_vertex_values(mesh)
 
-    vFile.close()
     rhoFile.close()
     stressFile.close()
 
@@ -68,19 +62,19 @@ def visualize(params, DIR=''):
     fig, axes = plt.subplots(2, 1, sharex=True, figsize=(8,8))
     axes[-1].set_xlabel(r'$x$')
     axes[0].set_ylabel(r'$\rho/\rho_0$')
-    axes[1].set_ylabel(r'$v/v_0$')
+    axes[1].set_ylabel('stress')
     axes[0].set_xlim(np.min(geometry), np.max(geometry))
     axes[0].set_ylim(0, np.max(rho)+0.05)
-    axes[1].set_ylim(np.min(v)-0.05, np.max(v)+0.05)
+    axes[1].set_ylim(np.min(stress)-0.05, np.max(stress)+0.05)
     
     rhoplot, = axes[0].plot(geometry[0], rho[0], 'g-', ms=3)
-    velplot, = axes[1].plot(geometry[0], v[0],   'r-', ms=3)
+    velplot, = axes[1].plot(geometry[0], stress[0],   'r-', ms=3)
     
     def update(value):
         ti = np.abs(times - value).argmin()
         rhoplot.set_ydata(rho[ti])
         rhoplot.set_xdata(geometry[ti])
-        velplot.set_ydata(v[ti])
+        velplot.set_ydata(stress[ti])
         velplot.set_xdata(geometry[ti])
         plt.draw()
         
@@ -91,23 +85,23 @@ def visualize(params, DIR=''):
     slider.drawon = False
     slider.on_changed(update)
 
-    # plt.show()  
-    print('Saving movie-...')
-    FPS = 100
-    movFile = os.path.join(DIR, '%s_fields.mov' % params['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 progressbar.progressbar(range(len(times))):
-        slider.set_val(times[tt])
-        fr = get_filename("%03d.png" % tt)
-        fig.savefig(fr, facecolor=fig.get_facecolor(), dpi=100)
-    os.system(command + " " + str(fps)
-                + " -i " + tmp_dir.name + os.sep
-                + "%03d.png " + options + " " + movFile)
-    tmp_dir.cleanup()
+    plt.show()  
+    # print('Saving movie-...')
+    # FPS = 100
+    # movFile = os.path.join(DIR, '%s_fields.mov' % params['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 progressbar.progressbar(range(len(times))):
+    #     slider.set_val(times[tt])
+    #     fr = get_filename("%03d.png" % tt)
+    #     fig.savefig(fr, facecolor=fig.get_facecolor(), dpi=100)
+    # os.system(command + " " + str(fps)
+    #             + " -i " + tmp_dir.name + os.sep
+    #             + "%03d.png " + options + " " + movFile)
+    # tmp_dir.cleanup()
 
     # plotting the stresses
     active_stress = -params['lamda'] * rho/(rho + params['saturation_rho'])
@@ -139,39 +133,39 @@ def visualize(params, DIR=''):
             fc='#999999')
     slider.drawon = False
     slider.on_changed(update_2)
-    # plt.show()
+    plt.show()
 
     # make movie
-    print('Saving movie-...')
-    FPS = 50
-    movFile = os.path.join(DIR, '%s_stresses.mov' % params['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 progressbar.progressbar(range(len(times))):
-        slider.set_val(times[tt])
-        fr = get_filename("%03d.png" % tt)
-        fig1.savefig(fr, facecolor=fig1.get_facecolor(), dpi=100)
-    os.system(command + " " + str(fps)
-                + " -i " + tmp_dir.name + os.sep
-                + "%03d.png " + options + " " + movFile)
-    tmp_dir.cleanup()
+    # print('Saving movie-...')
+    # FPS = 50
+    # movFile = os.path.join(DIR, '%s_stresses.mov' % params['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 progressbar.progressbar(range(len(times))):
+    #     slider.set_val(times[tt])
+    #     fr = get_filename("%03d.png" % tt)
+    #     fig1.savefig(fr, facecolor=fig1.get_facecolor(), dpi=100)
+    # os.system(command + " " + str(fps)
+    #             + " -i " + tmp_dir.name + os.sep
+    #             + "%03d.png " + options + " " + movFile)
+    # tmp_dir.cleanup()
 
     # L(t) vs t
-    length = geometry[:,-1,0]-geometry[:,0,0]
-    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)-0.05, np.max(length)+0.05)
-    ax2.plot(times, length)
-    plt.savefig("%s.png" %params['timestamp'])
-    # plt.show()     
-    np.save('%s_length.npy' %params['timestamp'], length)
-    print(length[-1])
-    print(length[-1]-length[-3])
+    # length = geometry[:,-1,0]-geometry[:,0,0]
+    # 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)-0.05, np.max(length)+0.05)
+    # ax2.plot(times, length)
+    # plt.savefig("%s.png" %params['timestamp'])
+    # # plt.show()     
+    # np.save('%s_length.npy' %params['timestamp'], length)
+    # print(length[-1])
+    # print(length[-1]-length[-3])
 
     # fig3, ax3 = plt.subplots(1, 1, figsize=(8, 8))
     # ax3.set_xlabel('$\log t$')