ditch density field altogether

u_v_rho/1D/check_response_curves_without_activity/parameters.json

@@ -12,6 +12,6 @@
     "timestep": 0.01,
     "turnover_rho": 1,
     "viscosity": 1,
-    "applied_boundary_stress": 2.0,
+    "applied_boundary_stress": 1.5,
     "zero_rho_boundary": true
 }
\ No newline at end of file

u_v_rho/1D/check_response_curves_without_activity/run_tissue.py

@@ -25,7 +25,6 @@ if args.jsonfile=='parameters.json':
     parametersFile = timestamp + '_parameters.json'
     initu = None
     initv = None
-    initrho = None
     initTime = 0.0
     mesh = None
 
@@ -74,7 +73,7 @@ else:
     parameters['maxtime'] = args.time
 
 tissue = OneDimTissue(parameters, mesh)
-tissue.solve(initu, initv, initrho, initTime, extend)
+tissue.solve(initu, initv, initTime, extend)
 
 if extend:
     parameters['maxtime'] = initTime + parameters['maxtime']

u_v_rho/1D/check_response_curves_without_activity/tissue.py

@@ -37,25 +37,14 @@ class OneDimTissue(object):
 
     def setup(self):
 
-        scalar_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1)
         vector_element = df.VectorElement('P', self.mesh.ufl_cell(), 1)
 
-        # u, v, rho
-        self.mixed_element = df.MixedElement([vector_element, vector_element, scalar_element])
+        # u, v
+        self.mixed_element = df.MixedElement([vector_element, vector_element])
 
         # define function space with this mixed element
         self.function_space = df.FunctionSpace(self.mesh, self.mixed_element)
         
-        # boundary condition on rho
-        if self.zero_rho_boundary:
-            self.rho_boundary = 0.0
-        else:
-            self.rho_boundary = self.average_rho
-
-        self.bc = df.DirichletBC(self.function_space.sub(2), 
-                                 df.Constant(self.rho_boundary), 
-                                 'on_boundary')
-        
         # define functions on function space
         self.function0 = df.Function(self.function_space)
         self.function  = df.Function(self.function_space)
@@ -100,30 +89,11 @@ class OneDimTissue(object):
         viscous_stress = self.viscosity * 2 * eps_v
         return (elastic_stress + viscous_stress)
 
-    def diffusion_reaction_rho(self, rho, trho):
-        return (self.diffusion_rho * df.inner(df.nabla_grad(rho), 
-                                              df.nabla_grad(trho)) 
-                - self.turnover_rho * df.inner(rho*(1 - rho/self.average_rho), 
-                                               trho)
-                )
-
     def setup_initial_conditions(self):
         # ic for u
         zero_vector = df.Constant((0.0,))
         u0 = df.interpolate(zero_vector, 
                             self.function_space.sub(0).collapse())
-        # ic for rho
-        if self.zero_rho_boundary:
-            base_rho = 0.0
-        else:
-            base_rho = self.average_rho
-        rho0 = df.interpolate(df.Expression(
-                            'base_rho + 0.1 * 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(2).collapse())
-
         
         # velocity from u0 and rho0
         VFS = self.function_space.sub(1).collapse()
@@ -137,12 +107,12 @@ class OneDimTissue(object):
                 - df.inner(tv0, self.imposed_boundary_stress_normal_component()) * self.ds(2)
         )
         df.solve(v0form == 0, v0)
-        df.assign(self.function0, [u0, v0, rho0])
+        df.assign(self.function0, [u0, v0])
 
     def setup_weak_forms(self):
-        u0, v0, rho0  = df.split(self.function0)
-        u,  v, rho    = df.split(self.function)
-        tu, tv, trho  = df.TestFunctions(self.function_space)
+        u0, v0  = df.split(self.function0)
+        u,  v    = df.split(self.function)
+        tu, tv  = df.TestFunctions(self.function_space)
                   
         uform = (df.inner((u - u0)/self.timestep, tu)
                 - df.inner(v0, tu))* df.dx
@@ -152,14 +122,10 @@ class OneDimTissue(object):
                            self.epsilon(tv))) * df.dx
                 + df.inner(tv, self.imposed_boundary_stress_normal_component()) * self.ds(1)
                 - df.inner(tv, self.imposed_boundary_stress_normal_component()) * self.ds(2))
-
-        rhoform = (df.inner((rho - rho0)/self.timestep, trho)
-                + self.advection(rho0, v0, trho)
-                + self.diffusion_reaction_rho(rho, trho))* df.dx
         
-        self.form = (uform + vform + rhoform)
+        self.form = (uform + vform)
 
-    def solve(self, initu=None, initv=None, initrho=None, initTime=0, extend=False):
+    def solve(self, initu=None, initv=None, initTime=0, extend=False):
 
         # create function space, functions, bc
         self.setup()
@@ -167,7 +133,6 @@ class OneDimTissue(object):
         # create files if they don't exist, open them if they do exist
         self.uFile   = df.XDMFFile('%s_displacement.xdmf' % self.timestamp)
         self.vFile   = df.XDMFFile('%s_velocity.xdmf' % self.timestamp)
-        self.rhoFile = df.XDMFFile('%s_density.xdmf' % self.timestamp)
         
         # time-variables
         self.time = initTime
@@ -179,10 +144,9 @@ class OneDimTissue(object):
             self.setup_initial_conditions()
 
             # save the ic
-            u0, v0, rho0 = self.function0.split(deepcopy=True)
+            u0, v0 = self.function0.split(deepcopy=True)
             self.uFile.write_checkpoint(u0, 'displacement', self.time)
             self.vFile.write_checkpoint(v0, 'velocity', self.time)
-            self.rhoFile.write_checkpoint(rho0, 'density', self.time)
 
             # move the mesh with this initial velocity
             dr0 = df.project(v0 * self.timestep, self.function_space.sub(0).collapse())
@@ -192,32 +156,26 @@ class OneDimTissue(object):
             #  assign fields to func0
             u0 = df.project(initu, self.function_space.sub(0).collapse())
             v0 = df.project(initv, self.function_space.sub(1).collapse())
-            rho0 = df.project(initrho, self.function_space.sub(2).collapse())
-            df.assign(self.function0, [u0, v0, rho0])
+            df.assign(self.function0, [u0, v0])
 
             # move the mesh with initial v read from file
             dr0 = df.project(v0 * self.timestep, self.function_space.sub(0).collapse())
             df.ALE.move(self.mesh, dr0)
 
-            self.setup()  # Builds space on moved, unrefined mesh
-            self.function0 = df.project(self.function0, self.function_space)
-
-
         self.setup_weak_forms()
 
         for steps in progressbar.ProgressBar()(range(1, maxsteps + 1)):
             self.time += self.timestep  
 
             # solve to get fields at next timestep
-            df.solve(self.form == 0, self.function, self.bc)
+            df.solve(self.form == 0, self.function)
 
             # save the fields
-            u, v, rho = self.function.split(deepcopy=True)
+            u, v = self.function.split(deepcopy=True)
             if steps % savesteps == 0:
 
                 self.uFile.write_checkpoint(u, 'displacement', self.time, append=True)
                 self.vFile.write_checkpoint(v, 'velocity', self.time, append=True)
-                self.rhoFile.write_checkpoint(rho, 'density', self.time, append=True)
 
             # assign the calculated function to a newly-defined function, to be used later
             old_function = df.Function(self.function_space)
@@ -241,5 +199,4 @@ class OneDimTissue(object):
             
         self.uFile.close()
         self.vFile.close()
-        self.rhoFile.close()
         
\ No newline at end of file

u_v_rho/1D/check_response_curves_without_activity/viz_tissue.py

@@ -16,7 +16,7 @@ def visualize(params, DIR=''):
     times = np.arange(savesteps + 1) * params['savetime']
     
     # Read mesh geometry and topology from h5 file
-    var = 'density'
+    var = 'displacement'
     h5 = h5py.File(os.path.join(DIR, '%s_%s.h5' % (params['timestamp'], var)), "r")
     
     geometry = []
@@ -32,11 +32,9 @@ def visualize(params, DIR=''):
     # Initialize field storage
     u = []
     v = []
-    rho = []
 
     uFile = df.XDMFFile(os.path.join(DIR, '%s_displacement.xdmf' % params['timestamp']))
     vFile = df.XDMFFile(os.path.join(DIR, '%s_velocity.xdmf' % params['timestamp']))
-    rhoFile = df.XDMFFile(os.path.join(DIR, '%s_density.xdmf' % params['timestamp']))
 
     for k in progressbar.ProgressBar()(range(len(times))):
         # Create and process the mesh
@@ -59,21 +57,18 @@ def visualize(params, DIR=''):
         # Read field values for the current mesh
         SFS = df.FunctionSpace(mesh, 'P', 1)
         VFS = df.VectorFunctionSpace(mesh, 'P', 1)
-        ui, vi, rhoi = df.Function(VFS), df.Function(VFS), df.Function(SFS)
+        ui, vi = df.Function(VFS), df.Function(VFS)
 
         uFile.read_checkpoint(ui, 'displacement', k)
         vFile.read_checkpoint(vi, 'velocity', k)
-        rhoFile.read_checkpoint(rhoi, 'density', k)
 
         u.append(ui.compute_vertex_values(mesh))
         v.append(vi.compute_vertex_values(mesh))
-        rho.append(rhoi.compute_vertex_values(mesh))
 
         del mesh
 
     uFile.close()
     vFile.close()
-    rhoFile.close()
 
     # Sort the geometry and fields
     for k in range(len(times)):
@@ -81,32 +76,26 @@ def visualize(params, DIR=''):
         geometry[k] = geometry[k][sorted_indices]
         u[k] = u[k][sorted_indices]
         v[k] = v[k][sorted_indices]
-        rho[k] = rho[k][sorted_indices]
 
     # minima and maxima of the fields
-    min_rho, min_u, min_v = map(lambda x: min(min(sublist) for sublist in x), [rho, u, v])
-    max_rho, max_u, max_v = map(lambda x: max(max(sublist) for sublist in x), [rho, u, v])
+    min_u, min_v = map(lambda x: min(min(sublist) for sublist in x), [u, v])
+    max_u, max_v = map(lambda x: max(max(sublist) for sublist in x), [u, v])
     min_geometry = min(min(sublist) for sublist in geometry)
     max_geometry = max(max(sublist) for sublist in geometry)
     # interactive plot
-    fig, axes = plt.subplots(3, 1, sharex=True, figsize=(8,8))
+    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/\ell \kappa$')
-    axes[2].set_ylabel(r'$u/\ell$')
+    axes[0].set_ylabel(r'$v/\ell \kappa$')
+    axes[1].set_ylabel(r'$u/\ell$')
     axes[0].set_xlim(min_geometry, max_geometry)
-    axes[0].set_ylim(min_rho-0.1, max_rho+0.1)
-    axes[1].set_ylim(min_v-0.1, max_v+0.1)
-    axes[2].set_ylim(min_u-0.1, max_u+0.1)
+    axes[0].set_ylim(min_v-0.1, max_v+0.1)
+    axes[1].set_ylim(min_u-0.1, max_u+0.1)
     
-    rhoplot, = axes[0].plot(geometry[0], rho[0], ms=3, color='g', lw = 2, marker = 'o')
-    velplot, = axes[1].plot(geometry[0], v[0], ms=3, color='r', lw = 2, marker = 'o')
-    uplot, = axes[2].plot(geometry[0], u[0], ms=3, color='b', lw = 2, marker = 'o')
+    velplot, = axes[0].plot(geometry[0], v[0], ms=3, color='r', lw = 2, marker = 'o')
+    uplot, = axes[1].plot(geometry[0], u[0], ms=3, color='b', lw = 2, marker = 'o')
     
     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_xdata(geometry[ti])
         uplot.set_ydata(u[ti])