Deleted growth_eqns_moving_mesh.py

growth_eqns_moving_mesh.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.f0 = df.Function(self.function_space)                 # f at t =0
-
-        df.plot(self.mesh)
-        plt.show()
-
-    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, trho):
-        return df.inner(self.turnover_rho*(rho - self.mean_density), trho)
-
-    def reaction_diffusion_c(self,c,tc):
-        return ( self.diffusion_c*df.inner(c.dx(0),tc.dx(0))
-                + self.reaction_c(c,tc)
-                )
-    def reaction_diffusion_rho(self,c,tc):
-        return ( self.diffusion_rho*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])
-
-    def setup_weak_forms(self):
-        u0, v0, 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)
-                + self.advection(c, v, tc)
-                + self.reaction_diffusion_c(c, tc)
-            )
-        rhoform = (df.inner((rho-rho0)/self.timestep,trho)
-                + self.advection(rho,v,trho)
-                + self.reaction_diffusion_rho(rho, 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.f0.assign(self.f)
-
-            df.ALE.move(self.mesh,u)
-        df.plot(self.mesh)
-        plt.show()
-        return (u_array, v_array, rho_array, c_array)
-
-
-
-if __name__ == '__main__':
-    import dolfin as df
-
-    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('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