given velocity, find whether rho=0 is a stable or unstable fixed point

basic_model_only_density/density_given_velocity.py

+import dolfin as df
+import numpy as np
+import progressbar
+import scipy as sc
+import os
+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])
+
+        # define the growth velocity field
+        self.sigma = df.Expression('x[0]*exp(-t)', t=0, degree=1)
+        # set up mesh
+        self.mesh = df.IntervalMesh(self.resolution, -self.system_size/2, self.system_size/2)
+        
+        # create mesh, function space, define function, test function
+        self.SFS = df.FunctionSpace(self.mesh, 'P', 1)
+        self.sfs1 = df.FunctionSpace(self.mesh, 'P', 1)
+        self.rho = df.Function(self.SFS)
+        self.rho0 = df.Function(self.SFS) 
+        trho = df.TestFunction(self.SFS)
+        self.bc = df.DirichletBC(self.SFS, df.Constant(0.0), 'on_boundary')
+        
+        self.velocity = df.project(self.sigma, self.sfs1)
+        self.form = (df.inner((self.rho - self.rho0)/self.timestep, trho)
+                   + df.inner((self.rho0 * self.velocity).dx(0), trho)
+                   + self.diffusion_rho * df.inner(self.rho.dx(0), trho.dx(0))
+                   - self.turnover_rho * df.inner(self.rho * (1 - self.rho/self.average_rho), trho)
+                    ) * df.dx
+
+
+    def solve(self):
+        savesteps = int(self.savetime/self.timestep)           
+        maxsteps = int(self.maxtime/self.timestep)
+        self.time = 0.0
+        rhoFile = df.XDMFFile('%s_density.xdmf' %params['timestamp'])   # create the file here
+
+        # initial condition
+
+        rho0 = df.Expression('cos(pi * x[0]/L) * cos(pi * x[0]/L)', pi = np.pi, L = self.system_size,  degree=1)
+        self.rho0.assign(df.project(rho0, self.SFS))
+
+        # save data
+        rhoFile.write_checkpoint(self.rho0, 'density', self.time)
+
+        # time stepping
+        for steps in progressbar.progressbar(range(1, maxsteps + 1)):
+            # get velocity
+            self.sigma.t = self.time
+            self.velocity.assign(df.project(self.sigma, self.sfs1))
+            # solve
+            df.solve(self.form == 0, self.rho, self.bc)
+            # update
+            self.rho0.assign(self.rho)
+            # save data
+            if steps % savesteps == 0:
+                rhoFile.write_checkpoint(self.rho0, 'density', self.time, append=True)
+            # move mesh
+            displacement = df.project(self.velocity * self.timestep, self.sfs1)
+            df.ALE.move(self.mesh, displacement)
+            self.time += self.timestep
+        rhoFile.close()
+    
+
+if __name__ == '__main__':
+    import json, datetime
+    assert os.path.isfile('parameters.json'), 'parameters.json file not found'  # assert that parameters.json is a valid file, otherwise 
+                                                                                # give an error message parameters.json file not found
+    
+    # load the parameters
+    with open('parameters.json') as jsonFile:                    
+        params = json.load(jsonFile)
+           
+    timestamp = datetime.datetime.now().strftime("%d%m%y-%H%M%S")
+    params['timestamp'] = timestamp
+        
+    gd = Growth(params)
+    gd.solve()
+    
+    with open(params['timestamp'] + '_parameters.json', "w") as fp:
+        json.dump(params, fp, indent=4) 
+    from viz_density_given_velocity import visualize
+    visualize(params, DIR='')  

basic_model_only_density/viz_density_given_velocity.py

+import dolfin as df
+import numpy as np
+import vedo as vd
+import os
+import h5py
+from matplotlib.widgets import Slider
+import progressbar
+import matplotlib.pyplot as plt
+from tempfile import TemporaryDirectory
+import numdifftools as nd
+
+def visualize(params, DIR=''):
+    savesteps = int(params['maxtime']/params['savetime'])
+    times = np.arange(savesteps+1) * params['savetime']    
+    # Read mesh geometry from h5 file
+    var = 'density'
+    h5 = h5py.File(os.path.join(DIR, 
+                                '%s_%s.h5' % (params['timestamp'], var)), "r")
+    # should be in the loop if remeshing
+    topology = np.array(h5['%s/%s_0/mesh/topology'%(var,var)])
+    geometry = []
+    for i in range(len(times)):
+        geometry.append(np.array(h5['%s/%s_%d/mesh/geometry'%(var,var,i)]))
+    h5.close()
+    geometry = np.array(geometry)
+    geometry, zeros= np.dsplit(geometry, 2)
+
+    mesh = df.IntervalMesh(params['resolution'], 
+                            - params['system_size']/2, 
+                              params['system_size']/2)
+    
+    # Read data
+    rho = np.zeros((len(times), mesh.num_vertices()))
+
+    rhoFile = df.XDMFFile(os.path.join(DIR, 
+                                    '%s_density.xdmf' % params['timestamp']))
+
+    # Reading data
+    print('Reading data...')
+    for steps in progressbar.progressbar(range(savesteps+1)):
+        mesh.coordinates()[:] = geometry[steps]
+        SFS = df.FunctionSpace(mesh, 'P', 1)
+        rhoi = df.Function(SFS)
+            
+        rhoFile.read_checkpoint(rhoi, 'density', steps)
+
+        rho[steps] = rhoi.compute_vertex_values(mesh)
+
+    rhoFile.close()
+
+    # interactive plot
+    plt.rcParams.update({'font.size': 18})
+    plt.yticks(fontsize=18)
+    plt.xticks(fontsize=18)
+    fig, axes = plt.subplots(1, 1, sharex=True, figsize=(8,8))
+    axes.set_xlabel(r'$x$')
+    axes.set_ylabel(r'$\rho/\rho_0$')
+    axes.set_xlim(np.min(geometry), np.max(geometry))
+    axes.set_ylim(0, np.max(rho)+0.05)
+    
+    rhoplot, = axes.plot(geometry[0], rho[0], 'g-', ms=3)
+    
+    def update(value):
+        ti = np.abs(times - value).argmin()
+        rhoplot.set_ydata(rho[ti])
+        rhoplot.set_xdata(geometry[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()  
+    # 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()
+
+if __name__ == '__main__':  
+    import argparse, json
+    parser = argparse.ArgumentParser()
+    parser.add_argument('-j','--jsonfile', help='parameters file', required=True)
+    parser.add_argument('--onscreen', action=argparse.BooleanOptionalAction)
+    args = parser.parse_args()
+
+    with open(args.jsonfile) as jsonFile:
+        params = json.load(jsonFile)
+
+    visualize(params, DIR=os.path.dirname(args.jsonfile))
\ No newline at end of file