recovered bois limit

fixed_boundaries.py

+import dolfin as df
+import numpy as np
+import progressbar
+
+df.set_log_level(df.LogLevel.ERROR)
+df.parameters['form_compiler']['optimize'] = True
+
+class FixedBoundaries(object):
+    def __init__(self, parameters):
+        # read in parameters
+        for key in parameters:
+            setattr(self, key, parameters[key])
+
+        # create mesh, mixed element and function space
+        self.mesh = df.IntervalMesh(self.resolution, 0.0, 2.0*np.pi*self.system_size_by_2PI)
+        scalar_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1)
+        mixed_element = df.MixedElement([scalar_element, scalar_element, scalar_element]) # u, v, rho
+
+        # define function space with this mixed element
+        self.function_space = df.FunctionSpace(self.mesh, mixed_element)
+
+        # dirichlet boundaries for u, v
+        bc1 = df.DirichletBC(self.function_space.sub(0), df.Constant(0.0), 'on_boundary')   # u,v at boundary = 0
+        bc2 = df.DirichletBC(self.function_space.sub(1), df.Constant(0.0), 'on_boundary')
+        self.bc = ([bc1,bc2])
+
+        # define the reqd functions on this function space
+        self.f  = df.Function(self.function_space)                # f at current time
+        self.f1 = df.Function(self.function_space)                # f at t = -1
+        self.f0 = df.Function(self.function_space)                # f at t =0
+
+    def advection(self, conc, vel, tconc):
+        return (df.inner((vel * conc).dx(0), tconc))
+
+    def active_stress(self, rho):
+        return self.lamda * rho/(rho + self.saturation_rho)
+    
+    def passive_stress(self, u, v):
+        return (self.elasticity * u.dx(0) + self.viscosity * v.dx(0))
+    
+    def stress(self, u, v, rho):
+        return (self.passive_stress(u, v) + self.active_stress(rho))
+        
+    def reaction_rho(self, rho, trho):
+        return self.turnover_rho * df.inner(rho - self.average_rho, trho)
+
+    def diffusion_reaction_rho(self, rho, trho):
+        return (self.diffusion_rho * df.inner(rho.dx(0), trho.dx(0)) + self.reaction_rho(rho, trho))
+
+    def setup_initial_conditions(self, u0, rho0):
+        u0   = df.interpolate(u0, self.function_space.sub(0).collapse())
+        rho0 = df.interpolate(rho0, self.function_space.sub(2).collapse())
+        VFS = self.function_space.sub(1).collapse()
+        v  = df.Function(VFS)
+        tv = df.TestFunction(VFS)
+        vform = (self.friction * df.inner(v, tv) 
+                    + df.inner(self.stress(u0, v, rho0), tv.dx(0))
+                ) * df.dx
+        bc = df.DirichletBC(VFS, df.Constant(0.0), 'on_boundary')
+        df.solve(vform == 0, v, bc)
+        df.assign(self.f0, [u0, v, rho0])
+        self.f1.assign(self.f0)
+
+    def setup_weak_forms(self):
+        u1, v1, rho1  = df.split(self.f1)
+        u0, v0, rho0  = df.split(self.f0)
+        u,  v, rho    = df.split(self.f)
+        tu, tv, trho = df.TestFunctions(self.function_space)
+
+        uform = (df.inner((u - u0)/self.timestep, tu)
+                - 9/16 * df.inner(v, tu)
+                - 3/8 * df.inner(v0, tu)
+                - 1/16 * df.inner(v1, tu)
+                )
+        
+        vform = (self.friction * df.inner(v, tv) 
+                + df.inner(self.stress(u, v, rho), tv.dx(0))
+                )
+
+        rhoform = (df.inner((rho - rho0)/self.timestep, trho)
+                + 3/2 * self.advection(rho0, v, trho)
+                - 1/2 * self.advection(rho1, v, trho)
+                + 9/16 * self.diffusion_reaction_rho(rho, trho)
+                + 3/8  * self.diffusion_reaction_rho(rho0, trho)
+                + 1/16 * self.diffusion_reaction_rho(rho1, trho)
+        )
+        self.form = (uform + vform + rhoform) * df.dx
+
+
+    def solve(self, u0, rho0):
+        times = np.arange(0.0, self.maxtime, self.timestep)
+        x = self.mesh.coordinates()
+        u_array = np.zeros((len(times), len(x)))
+        v_array = np.zeros_like(u_array)
+        rho_array = np.zeros_like(u_array)
+
+        self.setup_initial_conditions(u0, rho0)
+        self.setup_weak_forms()
+
+        u, v, rho = self.f0.split(deepcopy=True)
+        u_array[0] = u.compute_vertex_values(self.mesh)
+        v_array[0] = v.compute_vertex_values(self.mesh) 
+        rho_array[0] = rho.compute_vertex_values(self.mesh)
+    
+        for i in progressbar.progressbar(range(0, len(times))):
+            df.solve(self.form == 0, self.f, self.bc)
+            u, v, rho = self.f.split(deepcopy=True)
+            u_array[i] = u.compute_vertex_values(self.mesh)
+            v_array[i] = v.compute_vertex_values(self.mesh)
+            rho_array[i] = rho.compute_vertex_values(self.mesh)
+            self.f0.assign(self.f)
+            self.f1.assign(self.f0)
+        return (u_array, v_array, rho_array)
+
+
+if __name__ == '__main__':
+    import json
+    import matplotlib.pyplot as plt
+    from matplotlib.widgets import Slider
+    
+    with open('growth_parameters.json') as jsonFile:
+        params = json.load(jsonFile)
+
+    fb = FixedBoundaries(params)
+
+    u0 = df.Expression('sin(x[0])', degree=1)
+    rho0 = df.Expression('rho0 + 0.1 * (cos(x[0])+cos(x[0]/2))', rho0=fb.average_rho, degree=1) 
+
+    (u, v, rho) = fb.solve(u0, rho0)
+
+
+    x = fb.mesh.coordinates()
+    times = np.arange(0, fb.maxtime, fb.timestep)    
+    fig, (axu, axv, axrho) = plt.subplots(3,1, sharex=True, figsize=(8,6))
+    axrho.set_xlabel(r'$x$')
+    axu.set_ylabel(r'$u(x,t)$')
+    axv.set_ylabel(r'$v(x,t)$')
+    axrho.set_ylabel(r'$\rho(x,t)$')
+    axu.set_ylim(np.min(u), np.max(u))
+    axv.set_ylim(np.min(v), np.max(v))
+    axrho.set_ylim(np.min(rho), np.max(rho))
+    
+    uplot, = axu.plot(x, u[0])
+    vplot, = axv.plot(x, v[0])
+    rhoplot, = axrho.plot(x, rho[0])
+
+    def update(value):
+        ti = np.abs(times-value).argmin()
+        uplot.set_ydata(u[ti])
+        vplot.set_ydata(v[ti])
+        rhoplot.set_ydata(rho[ti])
+        plt.draw()
+        
+    sax = plt.axes([0.1, 0.92, 0.7, 0.02])
+    slider = Slider(sax, r'$t/\tau$', min(times), max(times),
+            valinit=min(times), valfmt='%3.1f',
+            fc='#999999')
+    slider.drawon = False
+    slider.on_changed(update)
+    plt.show()

growth_parameters.json

 {
-    "dimension" : 1,
-    "resolution" : 200,
+    "resolution" : 32,
     "system_size_by_2PI" : 1 ,
 
-    "elasticity" : 0.0,
+    "elasticity" : 1.0,
     "viscosity" : 1.0,
-    "lamda": 80.0,
+    "friction" : 1.0,
+    "lamda": 19.0,
 
     "diffusion_rho" : 1.0,
-    "turnover_rho" : 1.0,
+    "turnover_rho" : 0.0,
     "average_rho" : 1.0,
-    "saturation_rho" : 5.0,
+    "saturation_rho" : 1.0,
 
     "diffusion_c" : 1.0,
     "turnover_c" : 1.0,
     "average_c" : 1.0,
 
     "timestep" : 0.01,
-    "maxtime" : 100.0
+    "maxtime" : 10.0
 
 }