growth equations without f and F

growth_eqns.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 Growth(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,self.system_size)
+        conc_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1)
+        mixed_element = df.MixedElement([conc_element, conc_element, conc_element, conc_element]) # u, v, rho, c
+
+        # define function space with this mixed element
+        self.function_space = df.FunctionSpace(self.mesh,mixed_element)
+
+        # define the reqd functions on this function space
+        self.f = df.Function(self.function_space)                 # f at current time
+        self.fminus1 = df.Function(self.function_space)           # f at t = -1
+        self.f0 = df.Function(self.function_space)                 # f at t =0
+
+    def advection(self,c,v,tc):
+        return (df.inner((v*c).dx(0),tc))
+
+    def reaction_c(self,c,tc):
+        return df.inner(self.turnover_c*(c-self.mean_concentration),tc)
+
+    def reaction_rho(self,rho):
+        return self.turnover_rho*rho
+
+    def reaction_diffusion_c(self,c,tc):
+        return ( self.diffusion*df.inner(c.dx(0),tc.dx(0))
+                + self.reaction_c(c,tc)
+                )
+
+    def setup_initial_conditions(self, u0,rho0,c0):
+        
+        u0 = df.interpolate(u0, self.function_space.sub(0).collapse())
+
+        rho0 = df.interpolate(rho0, self.function_space.sub(2).collapse())
+
+        c0 = df.interpolate(c0, self.function_space.sub(3).collapse())
+
+        v0_function_space = df.FunctionSpace(self.mesh, 'P', 1)
+        v0  = df.Function(v0_function_space)
+        tv0 = df.TestFunction(v0_function_space)
+        v0form = (df.inner(v0,tv0)
+                + self.youngs_modulus*df.inner(u0.dx(0),tv0.dx(0))
+                - self.b*df.inner(rho0.dx(0),tv0)
+                ) * df.dx
+        df.solve(v0form==0, v0)
+        df.assign(self.f0, [u0, v0, rho0, c0])
+        self.fminus1.assign(self.f0)
+
+    def setup_weak_forms(self):
+        uminus1, _, rhominus1, cminus1  = df.split(self.fminus1)
+        u0, _ , rho0, c0 = df.split(self.f0)
+        u,  v, rho, c = df.split(self.f)
+        tu, tv, trho, tc = df.TestFunctions(self.function_space)
+
+        uform = (df.inner((u-u0)/self.timestep,tu) 
+                - df.inner(v,tu)
+        )
+
+        vform = (df.inner(v,tv)
+                + self.youngs_modulus*df.inner(u.dx(0),tv.dx(0))
+                - self.b*df.inner(rho.dx(0),tv)
+        )
+
+        cform = (df.inner((c-c0)/self.timestep, tc)
+                - 3/2 * self.advection(c0, v, tc)
+                + 1/2 * self.advection(cminus1, v, tc)
+                + 9/16 * self.reaction_diffusion_c(c, tc)
+                + 3/8 * self.reaction_diffusion_c(c0, tc)
+                + 1/16 * self.reaction_diffusion_c(cminus1, tc)
+            )
+
+        rhoform = (df.inner((rho-rho0)/self.timestep,trho)
+                - (3/2)*self.advection(rho0,v,trho)
+                + (1/2)*self.advection(rhominus1,v, trho)
+                + (9/16)*df.inner(self.reaction_rho(rho),trho)
+                + (3/8)*df.inner(self.reaction_rho(rho0),trho)
+                + (1/16)*df.inner(self.reaction_rho(rhominus1),trho)
+        )
+        self.form = (uform + vform + rhoform + cform) * df.dx
+
+    def solve(self,u0=None,rho0=None, c0=None):
+        times = np.arange(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)
+        c_array = np.zeros_like(u_array)
+
+        self.setup_initial_conditions(u0,rho0,c0)
+        self.setup_weak_forms()
+
+        u,v, rho, c = 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)
+        c_array[0] = c.compute_vertex_values(self.mesh)      
+
+        for i in progressbar.progressbar(range(0, len(times))):
+            df.solve(self.form == 0,self.f)
+            u, v, rho, c = 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)
+            c_array[i] = c.compute_vertex_values(self.mesh)
+            self.fminus1.assign(self.f0)
+            self.f0.assign(self.f)
+
+        return (u_array, v_array, rho_array, c_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)
+    assert(params['dimension']==1)
+
+    g = Growth(params)
+    u0 = df.Expression('0.5*cos(x[0]/2)', degree=1)  
+    rho0 = df.Expression('1+0.2*cos(x[0])', degree=1)
+    c0 = df.Expression('cos(2*x[0])+cos(x[0])', degree=1)
+    (u, v, rho, c) = g.solve(u0,rho0,c0)
+
+    x = g.mesh.coordinates()
+    times = np.arange(0,g.maxtime,g.timestep)    
+    fig, (axu, axv, axrho, axc) = plt.subplots(4,1, sharex=True, figsize=(8,6))
+    axu.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)$')
+    axc.set_ylabel(r'$c(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))
+    axc.set_ylim(np.min(c),np.max(c))
+    
+    uplot, = axu.plot(x, u[0])
+    vplot, = axv.plot(x, v[0])
+    rhoplot, = axrho.plot(x, rho[0])
+    cplot, = axc.plot(x, c[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])
+        cplot.set_ydata(c[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()
\ No newline at end of file