patterns on fixed boundaries

fixed_domain/fixed_bdry_fft.py

+import numpy as np
+from scipy.integrate import odeint
+import matplotlib.pyplot as plt
+from matplotlib.widgets import Slider
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+import time
+import progressbar 
+
+#parameters
+Lx = 2*np.pi
+Nx = 200
+T = 100
+dt = 0.01
+Nt = int(T/dt)
+times = np.linspace(0,T,Nt)
+
+K = 1
+eta = 1
+lamda = 5.0
+krho = 1
+rho0 = 1
+Drho = 1
+
+x = np.linspace(0, Lx, Nx)
+q = np.fft.fftshift(np.fft.fftfreq(len(x), d=Lx/(2*np.pi*Nx)))
+
+def reaction_rho(rho, rho0, krho):
+    return (-krho * (rho - rho0))
+
+def time_derivative(c, t):
+    # split
+    u, rho = np.split(c, 2)
+
+    # compute FFT
+    uq = np.fft.fftshift(np.fft.fft(u))
+    rhoq = np.fft.fftshift(np.fft.fft(rho - rho0))
+    vq = (-K * q**2 * uq + 1j * q * lamda * rhoq)/(1 + eta * q**2)
+
+    # RHS in Fourier-space
+    RHS_u_q = vq
+    RHS_rho_q = -rho0 * 1j * q * vq - Drho * q**2 * rhoq
+
+    # RHS in real space
+    RHS_u = np.real(np.fft.ifft(np.fft.ifftshift(RHS_u_q)))
+    RHS_rho = np.real(np.fft.ifft(np.fft.ifftshift(RHS_rho_q))) + reaction_rho(rho, rho0, krho)
+
+    return np.concatenate([RHS_u, RHS_rho])
+
+# initial conditions
+u00 = np.sin(x)
+# u00 = u00 * np.ones_like(x) + 0.01 * np.random.randn(x.shape[0])
+rho00 = np.cos(x) + np.cos(x/2)
+# rho00 = rho00 * np.ones_like(x) + 0.01 * np.random.randn(x.shape[0])
+
+c0 = np.concatenate([u00, rho00])
+
+#integrate in time
+c = odeint(time_derivative, c0, times)
+
+#split and reshape solution arrays
+u, rho = np.split(c, 2, axis=1)
+
+#plotting
+fig, (axu, axrho) = plt.subplots(2,1, sharex=True, figsize=(8,6))
+axu.set_xlabel(r'$x$')
+axu.set_ylabel(r'$u(x,t)$')
+axrho.set_ylabel(r'$\rho(x,t)$')
+axu.set_ylim(np.min(u), np.max(u))
+axrho.set_ylim(np.min(rho), np.max(rho))
+
+uplot, = axu.plot(x, u[0])
+rhoplot, = axrho.plot(x, rho[0])
+
+def update(value):
+    ti = np.abs(times-value).argmin()
+    uplot.set_ydata(u[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()
+

fixed_domain/fixed_boundaries.py

+import dolfin as df
+import numpy as np
+import progressbar
+import os
+
+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
+        L = 2.0*np.pi*self.system_size_by_2PI
+        if self.dimension==1:
+            self.mesh = df.IntervalMesh(self.resolution, 0.0, L)
+            assert self.bulk_elasticity == self.shear_elasticity
+            assert self.bulk_viscosity == self.shear_viscosity
+        elif self.dimension==2:
+            if hasattr(self, 'meshfile'):
+                print('Using %s' % str(self.meshfile))
+                self.mesh = df.Mesh(str(self.meshfile))
+            else:                
+                self.mesh = df.RectangleMesh(df.Point(0,0), df.Point(L,L), 
+                                            self.resolution, self.resolution)
+
+        scalar_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1)
+        vector_element = df.VectorElement('P', self.mesh.ufl_cell(), 1)
+        if self.morphogen:
+            # u, v, rho, c
+            mixed_element = df.MixedElement([vector_element, vector_element, 
+                                             scalar_element, scalar_element]) 
+            self.savedata = self.save_data_uvrhoc
+        else:
+            # u, v, rho
+            mixed_element = df.MixedElement([vector_element, vector_element, scalar_element])
+            self.savedata = self.save_data_uvrho
+
+        # define function space with this mixed element
+        self.function_space = df.FunctionSpace(self.mesh, mixed_element)
+
+        # dirichlet boundaries for u, v
+        # u,v at boundary = 0
+        if self.dimension==1:
+            self.zero_vector = df.Constant((0.0,))
+        elif self.dimension==2:
+            self.zero_vector = df.Constant((0.0, 0.0))
+        bc1 = df.DirichletBC(self.function_space.sub(0), self.zero_vector, 'on_boundary')
+        bc2 = df.DirichletBC(self.function_space.sub(1), self.zero_vector, '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(df.div(vel * conc), tconc))
+
+    def active_stress(self, rho, c):
+        return self.lamda * rho/(rho + self.saturation_rho) * c * df.Identity(self.dimension)
+    
+    def epsilon(self, v):
+        return df.sym(df.nabla_grad(v))
+    
+    def passive_stress(self, u, v):
+        eps_u, eps_v = self.epsilon(u), self.epsilon(v)
+        
+        elastic_stress = (self.shear_elasticity * eps_u + 
+                    (self.bulk_elasticity - self.shear_elasticity/self.dimension) * 
+                        df.tr(eps_u) * df.Identity(self.dimension))
+        
+        viscous_stress = (self.shear_viscosity * eps_v + 
+                    (self.bulk_viscosity - self.shear_viscosity/self.dimension) * 
+                        df.tr(eps_v) * df.Identity(self.dimension))
+        
+        return (elastic_stress + viscous_stress)
+    
+    def stress(self, u, v, rho, c):
+        return (self.passive_stress(u, v) + self.active_stress(rho, c))
+        
+    def reaction_rho(self, rho, trho):
+        return self.turnover_rho * df.inner(rho - self.average_rho, trho)
+
+    def reaction_c(self, c, tc):
+        return self.turnover_c * df.inner(c - self.average_c, tc)
+    
+    def diffusion_reaction_rho(self, rho, trho):
+        return (self.diffusion_rho * df.inner(df.nabla_grad(rho), df.nabla_grad(trho)) 
+                + self.reaction_rho(rho, trho))
+
+    def diffusion_reaction_c(self, c, tc):
+        return (self.diffusion_c * df.inner(df.nabla_grad(c), df.nabla_grad(tc)) 
+                    + self.reaction_c(c, tc))
+
+    def setup_initial_conditions(self, ui=None, rhoi=None, ci=None):
+        if (ui is not None) and (rhoi is not None):
+            u0   = df.interpolate(ui, self.function_space.sub(0).collapse())
+            rho0 = df.interpolate(rhoi, self.function_space.sub(2).collapse())
+            if self.morphogen:
+                c0 = df.interpolate(ci, self.function_space.sub(3).collapse())
+            else:
+                c0 = df.Constant(1.0)            
+        else:
+            u0   = df.interpolate(self.zero_vector, self.function_space.sub(0).collapse())
+            rho0 = df.interpolate(df.Constant(self.average_rho), self.function_space.sub(2).collapse())
+            # add noise
+            noise_u = self.noise_level * (2*np.random.random(u0.vector().size())-1)
+            u0.vector()[:] = u0.vector()[:] + noise_u            
+            noise_rho = self.noise_level * (2*np.random.random(rho0.vector().size())-1)
+            rho0.vector()[:] = rho0.vector()[:] + noise_rho
+            if self.morphogen:
+                c0 = df.interpolate(df.Constant(self.average_c), self.function_space.sub(3).collapse())
+                # add noise
+                noise_c = self.noise_level * (2*np.random.random(c0.vector().size())-1)
+                c0.vector()[:] = c0.vector()[:] + noise_c
+            else:
+                c0 = df.Constant(1.0)
+        
+        VFS = self.function_space.sub(1).collapse()
+        v0  = df.Function(VFS)
+        tv = df.TestFunction(VFS)
+        vform = (self.friction * df.inner(v0, tv) 
+                    + df.inner(self.stress(u0, v0, rho0, c0), self.epsilon(tv))
+                ) * df.dx
+        bc = df.DirichletBC(VFS, self.zero_vector, 'on_boundary')
+        df.solve(vform == 0, v0, bc)
+        if self.morphogen:
+            df.assign(self.f0, [u0, v0, rho0, c0])
+        else:
+            df.assign(self.f0, [u0, v0, rho0])
+        self.f1.assign(self.f0)
+
+    def setup_weak_forms(self):
+        if self.morphogen:
+            u1, v1, rho1, c1  = df.split(self.f1)
+            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)
+        else:
+            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)
+            c = df.Constant(1.0)
+                  
+        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, c), self.epsilon(tv))
+                )
+
+        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)
+                )
+        
+        if self.morphogen:
+            cform = (df.inner((c - c0)/self.timestep, tc)
+                + 3/2 * self.advection(c0, v, tc)
+                - 1/2 * self.advection(c1, v, tc)
+                + 9/16 * self.diffusion_reaction_c(c,  tc)
+                + 3/8  * self.diffusion_reaction_c(c0, tc)
+                + 1/16 * self.diffusion_reaction_c(c1, tc)
+                )
+            self.form = (uform + vform + rhoform + cform) * df.dx
+        else:
+            self.form = (uform + vform + rhoform) * df.dx
+            
+    def save_data_uvrho(self):
+        u, v, rho = self.f0.split(deepcopy=True)
+        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)
+        
+    def save_data_uvrhoc(self):
+        u, v, rho, c  = self.f0.split(deepcopy=True)
+        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)
+        self.cFile.write_checkpoint(c, 'concentration', self.time, append=True)
+        
+    def solve(self, extend=False, DIR=''):
+        # for saving data
+        self.uFile   = df.XDMFFile(os.path.join(DIR, '%s_displacement.xdmf' % self.timestamp))
+        self.vFile   = df.XDMFFile(os.path.join(DIR, '%s_velocity.xdmf' % self.timestamp))
+        self.rhoFile = df.XDMFFile(os.path.join(DIR, '%s_density.xdmf' % self.timestamp))
+        if self.morphogen:
+            self.cFile = df.XDMFFile(os.path.join(DIR, '%s_concentration.xdmf' % self.timestamp))
+        
+        if extend:
+            SFS = df.FunctionSpace(self.mesh, 'P', 1)
+            VFS = df.VectorFunctionSpace(self.mesh, 'P', 1)
+            ui, vi, rhoi, ci = df.Function(VFS), df.Function(VFS), df.Function(SFS), df.Function(SFS)
+            self.uFile.read_checkpoint(ui, 'displacement', -1)
+            self.vFile.read_checkpoint(vi, 'velocity', -1)
+            self.rhoFile.read_checkpoint(rhoi, 'density', -1)
+            if self.morphogen:
+                self.cFile.read_checkpoint(ci, 'concentration', -1)
+            else:
+                ci.interpolate(df.Constant(1.0))
+        else:
+            ui, vi, rhoi, ci = None, None, None, None     
+            
+        self.setup_initial_conditions(ui, rhoi, ci)
+        self.setup_weak_forms()
+        
+        # time-variables
+        self.time = 0.0
+        savesteps = int(self.savetime/self.timestep)
+        maxsteps = int(self.maxtime/self.timestep)
+        
+        if not extend:
+            self.savedata()
+        for steps in progressbar.progressbar(range(1, maxsteps + 1)):
+            df.solve(self.form == 0, self.f, self.bc)
+            self.time += self.timestep
+            df.assign(self.f1, self.f0)
+            df.assign(self.f0, self.f)
+            if steps % savesteps == 0:
+                self.savedata()
+
+        self.uFile.close()
+        self.vFile.close()
+        self.rhoFile.close()
+        if self.morphogen:
+            self.cFile.close()
\ No newline at end of file

fixed_domain/lsa/no_c_lsa.py

+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.widgets import Slider
+
+# real part of largest eigenvalue of STABILITY MATRIX
+def largest_real_eigval(q, gamma = 1.0, K=1.0, lamda=1.5, eta = 1, k_rho=10.0, rho_0=1, D_rho=1, rho_s=1.0):
+    A = np.asmatrix([[-K*q**2/(gamma+eta*q**2),1j*lamda*q*rho_s/((gamma+eta*q**2)*(rho_0 + rho_s)**2)], [1j*K*rho_0*q**3/(gamma+eta*q**2),-D_rho*q**2 -k_rho + rho_0*lamda*q**2*rho_s/((gamma+eta*q**2)*(rho_0 + rho_s)**2)]]) 
+    lamda = np.real(np.linalg.eigvals(A))
+    return lamda.max()
+
+kyu = np.linspace(0, 2*np.pi,1000)
+
+fig, ax = plt.subplots(1, figsize=(8,6))
+fig.subplots_adjust(left=0.15, bottom=0.3, right=0.98, top=0.95,
+        wspace=0.1, hspace=0.2)
+ax.set_xlabel(r'$q$')
+ax.set_ylabel(r'$Re[\, \lambda(q) \, ]_{\rm max}$')
+ax.axhline(y=0, color='black')
+
+gamma0 = 1.0
+D_rho_0 = 1.0
+lamda0 = 21.0
+K0 = 0.1
+k_rho0 = 0.1
+eta0 = 1.0
+rho_00 = 1.0
+rho_s0 = 1.0
+
+lamda = np.array([largest_real_eigval(q, gamma = gamma0, K=K0, lamda=lamda0, eta=eta0, k_rho=k_rho0, rho_0=rho_00, D_rho=D_rho_0, rho_s = rho_s0) for q in kyu])
+lambda_plot, = ax.plot(kyu, lamda)
+lamda_min, lamda_max = min(0, lamda.min()), max(0, lamda.max())
+ax.set_ylim(lamda_min, lamda_max)
+
+# axes for controlling parameter sliders
+ax_K  = plt.axes([0.1, 0.15, 0.2, 0.02])
+ax_lamda  = plt.axes([0.1, 0.20, 0.2, 0.02])
+ax_eta = plt.axes([0.4, 0.15, 0.2, 0.02])
+ax_k_rho = plt.axes([0.4, 0.20, 0.2, 0.02])
+ax_gamma  = plt.axes([0.4, 0.10, 0.2, 0.02])
+ax_rho_0  = plt.axes([0.7, 0.15, 0.2, 0.02])
+ax_D_rho  = plt.axes([0.7, 0.20, 0.2, 0.02])
+ax_rho_s  = plt.axes([0.7, 0.10, 0.2, 0.02])
+
+# sliders for controlling parameters
+s_K  = Slider(ax_K, r'$K$', valmin=0.0, valmax=50.0, valinit=K0, valstep=0.001)
+s_lamda  = Slider(ax_lamda, r'$\lambda$', valmin=0.0, valmax=100.0, valinit=lamda0, valstep=0.001)
+s_k_rho  = Slider(ax_k_rho, r'$k_{\rho}$', valmin=0.0, valmax=50.0, valinit=k_rho0, valstep=0.001)
+s_rho_0  = Slider(ax_rho_0, r'$\rho_0$', valmin=0.0, valmax=50.0, valinit=rho_00, valstep=0.001)
+s_eta  = Slider(ax_eta, r'$\eta$', valmin=0.0, valmax=50.0, valinit=eta0, valstep=0.001)
+s_D_rho  = Slider(ax_D_rho, r'$D_{\rho}$', valmin=0.0, valmax=50.0, valinit=D_rho_0, valstep=0.001)
+s_rho_s  = Slider(ax_rho_s, r'$\rho_s$', valmin=0.0, valmax=50.0, valinit=rho_s0, valstep=0.001)
+s_gamma  = Slider(ax_gamma, r'$\gamma$', valmin=0.0, valmax=50.0, valinit=gamma0, valstep=0.001)
+
+
+# slider update function
+def update(val):
+    lamda = np.array([largest_real_eigval(q, gamma = s_gamma.val, K=s_K.val, lamda=s_lamda.val, eta=s_eta.val, k_rho=s_eta.val, rho_0=s_rho_0.val, D_rho=s_D_rho.val, rho_s = s_rho_s.val)
+                            for q in kyu])
+    lambda_plot.set_ydata(lamda)
+    lamda_min, lamda_max = min(0, lamda.min()), max(0, lamda.max())
+    ax.set_ylim(lamda_min, lamda_max)
+    fig.canvas.draw()
+
+s_K.on_changed(update)
+s_lamda.on_changed(update)
+s_eta.on_changed(update)
+s_k_rho.on_changed(update)
+s_rho_0.on_changed(update)
+s_D_rho.on_changed(update)
+s_rho_s.on_changed(update)
+s_gamma.on_changed(update)
+
+plt.show()

fixed_domain/lsa/with_c_lsa.py

+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.widgets import Slider
+
+# real part of largest eigenvalue of STABILITY MATRIX
+def largest_real_eigval(q, K=1.0, lamda=1.5, eta = 1, krho=10.0, rho0 = 1, c0 = 1.0, kc = 1.0, Drho = 1.0, Dc = 1.0, gamma = 1.0):
+    A = np.asmatrix([[-K * q**2/(gamma + eta * q**2), 1j * c0 * q * lamda/(gamma + eta * q**2), 1j * rho0 * q * lamda/(gamma + eta * q**2)],
+                    [1j * q**3 * rho0 * K/(gamma + eta * q**2), -Drho * q**2 -krho + q**2 * c0 * rho0 * lamda/(gamma + eta * q**2), q**2 * rho0**2 * lamda/(gamma + eta * q**2)],
+                    [1j * q**3 * c0 * K/(gamma + eta * q**2), q**2 * c0**2 * lamda/(gamma + eta * q**2), -kc + q**2 * c0 * rho0 * lamda/(gamma + eta * q**2) - q**2 * Dc]]) 
+    lamda = np.real(np.linalg.eigvals(A))
+    return lamda.max()
+
+kyu = np.linspace(0,100,1000)
+
+fig, ax = plt.subplots(1, figsize=(8,6))
+fig.subplots_adjust(left=0.15, bottom=0.3, right=0.98, top=0.95,
+        wspace=0.1, hspace=0.2)
+ax.set_xlabel(r'$q$')
+ax.set_ylabel(r'$Re[\, \lambda(q) \, ]_{\rm max}$')
+ax.axhline(y=0, color='black')
+
+lamda0 = 2.6
+K0 = 1.0
+krho0 = 1.0
+eta0 = 1.0
+rho00 = 1.0
+Dc0 = 1.0
+Drho0 = 1.0
+kc0 = 1.0
+c00 = 1.0
+gamma0 = 1.0
+
+lamda = np.array([largest_real_eigval(q, K = K0, lamda = lamda0, eta = eta0, krho = krho0, rho0 = rho00, c0 = c00, kc = kc0, Drho = Drho0, Dc = Dc0, gamma = gamma0) for q in kyu])
+lambda_plot, = ax.plot(kyu, lamda)
+lamda_min, lamda_max = min(0, lamda.min()), max(0, lamda.max())
+ax.set_ylim(lamda_min, lamda_max)
+
+# axes for controlling parameter sliders
+ax_K  = plt.axes([0.1, 0.15, 0.2, 0.02])
+ax_lamda  = plt.axes([0.1, 0.20, 0.2, 0.02])
+ax_eta = plt.axes([0.4, 0.15, 0.2, 0.02])
+ax_krho = plt.axes([0.4, 0.20, 0.2, 0.02])
+ax_rho0  = plt.axes([0.7, 0.15, 0.2, 0.02])
+ax_Dc  = plt.axes([0.7, 0.20, 0.2, 0.02])
+ax_kc = plt.axes([0.1, 0.10, 0.2, 0.02])
+ax_c0 = plt.axes([0.4, 0.10, 0.2, 0.02])
+ax_Drho = plt.axes([0.7, 0.10, 0.2, 0.02])
+ax_gamma = plt.axes([0.7, 0.05, 0.2, 0.02])
+
+# sliders for controlling parameters
+s_K  = Slider(ax_K, r'$K$', valmin = 0.0, valmax = 10.0, valinit = K0, valstep = 0.001)
+s_lamda  = Slider(ax_lamda, r'$\lambda$', valmin = 0.0, valmax = 20.0, valinit = lamda0, valstep = 0.001)
+s_eta  = Slider(ax_eta, r'$\eta$', valmin = 0.0, valmax = 10.0, valinit = eta0, valstep = 0.001)
+s_krho  = Slider(ax_krho, r'$k_{\rho}$', valmin = 0.0, valmax = 10.0, valinit = krho0, valstep = 0.001)
+s_rho0  = Slider(ax_rho0, r'$\rho_0$', valmin = 0.0, valmax = 10.0, valinit = rho00, valstep = 0.001)
+s_Dc  = Slider(ax_Dc, r'$D_c$', valmin = 0.0, valmax = 10.0, valinit = Dc0, valstep = 0.001)
+s_kc  = Slider(ax_kc, r'$k_c$', valmin = 0.0, valmax = 10.0, valinit = kc0, valstep = 0.001)
+s_c0  = Slider(ax_c0, r'$c_0$', valmin = 0.0, valmax = 10.0, valinit = c00, valstep = 0.001)
+s_Drho  = Slider(ax_Drho, r'$D_\rho{}$', valmin = 0.0, valmax = 10.0, valinit = Drho0, valstep = 0.001)
+s_gamma  = Slider(ax_gamma, r'$\gamma$', valmin = 0.0, valmax = 10.0, valinit = gamma0, valstep = 0.001)
+
+# slider update function
+def update(val):
+    lamda = np.array([largest_real_eigval(q, K=s_K.val, lamda = s_lamda.val, eta = s_eta.val, krho = s_krho.val, rho0 = s_rho0.val, c0 = s_c0.val, kc = s_kc.val, Drho = s_Drho.val, Dc = s_Dc.val, gamma = s_gamma.val)
+                            for q in kyu])
+    lambda_plot.set_ydata(lamda)
+    lamda_min, lamda_max = min(0, lamda.min()), max(0, lamda.max())
+    ax.set_ylim(lamda_min, lamda_max)
+    fig.canvas.draw()
+
+s_K.on_changed(update)
+s_lamda.on_changed(update)
+s_eta.on_changed(update)
+s_krho.on_changed(update)
+s_rho0.on_changed(update)
+s_Dc.on_changed(update)
+s_Drho.on_changed(update)
+s_kc.on_changed(update)
+s_c0.on_changed(update)
+s_gamma.on_changed(update)
+
+plt.show()

fixed_domain/run_fixed_boundaries.py

+import datetime
+import json
+import os
+import argparse
+
+from fixed_boundaries import FixedBoundaries
+from viz_fixed_boundaries import visualize
+
+parser = argparse.ArgumentParser()
+parser.add_argument('-j','--jsonfile', help='json file', default='fixed_boundaries.json')
+parser.add_argument('-t','--time', help='time to extend', type=float, default=50)
+parser.add_argument('-o','--outputdir', help='output directory', type=str, default='')
+args = parser.parse_args()
+
+assert os.path.isfile(args.jsonfile), '%s file not found' % args.jsonfile
+with open(args.jsonfile) as jsonFile:
+    params = json.load(jsonFile)
+
+u0, rho0, c0 = None, None, None
+if args.jsonfile=='fixed_boundaries.json':
+    print('Fresh run')
+    extend = False
+    DIR = args.outputdir
+    timestamp = datetime.datetime.now().strftime("%d%m%y-%H%M%S")
+    # current timestamp is the jobID
+    params['timestamp'] = timestamp
+
+else:
+    print('Extending %s' % params['timestamp'])
+    extend = True
+    DIR = os.path.dirname(args.jsonfile)
+    oldMaxTime = params['maxtime']
+    params['maxtime'] = args.time
+
+# solve
+fb = FixedBoundaries(params)
+fb.solve(extend=extend, DIR=DIR)
+
+if extend:
+    params['maxtime'] = oldMaxTime + params['maxtime']
+        
+with open(os.path.join(DIR, 
+            params['timestamp'] + '_fixed_boundaries.json'), "w") as fp:
+    json.dump(params, fp, indent=4)
+
+visualize(params, DIR, interactive=True)
\ No newline at end of file

fixed_domain/scan_fixed_boundaries.py

+import datetime, time, json, os, itertools
+import dolfin as df
+import numpy as np
+import argparse
+from fixed_boundaries import FixedBoundaries
+from viz_fixed_boundaries import visualize
+
+parser = argparse.ArgumentParser()
+parser.add_argument('-o','--outputdir', help = 'output directory', type=str, default='')
+args = parser.parse_args()
+df.set_log_level(df.LogLevel.ERROR)
+df.parameters['form_compiler']['optimize'] = True
+
+# base parameters
+parameters = {
+    "morphogen" : False,
+
+    "resolution" : 32,
+    "system_size_by_2PI" : 1 ,
+
+    "elasticity" : 0.0,
+    "viscosity" : 1.0,
+    "friction" : 1.0,
+    "lamda": 5.5,
+
+    "diffusion_rho" : 1.0,
+    "turnover_rho" : 0.0,
+    "average_rho" : 1.0,
+    "saturation_rho" : 1.0,
+
+    "diffusion_c" : 0.0,
+    "turnover_c" : 0.0,
+    "average_c" : 0.0,
+
+    "noise_level": 0.1,
+
+    "timestep" : 0.01,
+    "savetime": 0.5,
+    "maxtime" : 200.0
+
+}
+
+# parameters to scan
+lamda_list = np.arange(1.0, 61.0, 20.0)
+turnover_rho_list = np.arange(0, 1.0, 0.1)
+elasticity_list = np.arange(0, 1.0, 0.1)
+
+plist = list(itertools.product(lamda_list, turnover_rho_list, elasticity_list))
+
+# scan parameters
+for i in range(len(plist)):
+    lamda, turnover_rho, elasticity = plist[i]
+
+    print('-'*50)
+    print('%d/%d' % (i, len(plist)))
+    print('lamda = %4.3f, turnover_rho = %4.3f, elasticity = %4.3f' % (lamda, turnover_rho, elasticity))
+    print('-'*50)
+
+    # timestamp as job-ID
+    timestamp = datetime.datetime.now().strftime("%d%m%y-%H%M%S")
+
+    params = parameters.copy()
+    params['timestamp'] = timestamp
+    params['turnover_rho'] = turnover_rho
+    params['elasticity'] = elasticity
+    params['lamda'] = lamda
+
+
+    try:
+        fb = FixedBoundaries(params)
+        DIR = args.outputdir
+
+        fb.solve(extend=False, DIR=DIR)
+        with open(os.path.join(DIR, params['timestamp']+'_fixed_boundaries.json'), "w") as fp:
+            json.dump(params, fp, indent=4)
+            
+        visualize(params, DIR=DIR, interactive=False)
+    except:
+        print('*** FAILED ****')
+        os.system("rm -f %s*" % timestamp)
+
+    time.sleep(2)

fixed_domain/viz_fixed_boundaries.py

+import numpy as np
+import json
+import matplotlib.pyplot as plt
+from matplotlib.widgets import Slider
+import argparse
+import progressbar
+import os
+import dolfin as df
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+import matplotlib.tri as tri
+from tempfile import TemporaryDirectory
+
+def mesh2triang(mesh):
+    xy = mesh.coordinates()
+    return tri.Triangulation(xy[:, 0], xy[:, 1], mesh.cells())
+
+def read_data(params, DIR):
+    print('Reading data from %s/%s...' % (DIR, params['timestamp']))
+    savesteps = int(params['maxtime']/params['savetime'])
+    times = np.arange(savesteps+1) * params['savetime']
+
+    # create mesh, mixed element and function space
+    L = 2.0*np.pi*params['system_size_by_2PI']
+    d = params['dimension']
+    if d==1:
+        mesh = df.IntervalMesh(params['resolution'], 0.0, L)
+    elif d==2:
+        if 'meshfile' in params.keys():
+            print('Using %s' % str(params['meshfile']))
+            mesh = df.Mesh(str(params['meshfile']))
+        else:
+            mesh = df.RectangleMesh(df.Point(0,0), df.Point(L,L), 
+                                    params['resolution'], params['resolution'])
+                
+    SFS = df.FunctionSpace(mesh, 'P', 1)
+    VFS = df.VectorFunctionSpace(mesh, 'P', 1)
+    ui, vi, rhoi = df.Function(VFS), df.Function(VFS), df.Function(SFS)
+    
+    u = np.zeros((len(times), mesh.num_vertices(), d))
+    v = np.zeros_like(u)
+    rho = np.zeros((len(times), mesh.num_vertices()))
+    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']))
+    if params['morphogen']:
+        ci = df.Function(SFS)
+        c = np.zeros_like(rho)
+        cFile = df.XDMFFile(os.path.join(DIR, '%s_concentration.xdmf' % params['timestamp']))
+    for steps in progressbar.progressbar(range(savesteps+1)):
+        uFile.read_checkpoint(ui, 'displacement', steps)
+        vFile.read_checkpoint(vi, 'velocity', steps)
+        rhoFile.read_checkpoint(rhoi, 'density', steps)
+        
+        u_vec = ui.compute_vertex_values(mesh)
+        u[steps] = u_vec.reshape(d, int(u_vec.shape[0]/d)).T
+
+        v_vec = vi.compute_vertex_values(mesh)
+        v[steps] = v_vec.reshape(d, int(v_vec.shape[0]/d)).T
+
+        rho[steps] = rhoi.compute_vertex_values(mesh)    
+        if params['morphogen']:
+            cFile.read_checkpoint(ci, 'concentration', steps)
+            c[steps] = ci.compute_vertex_values(mesh)
+    uFile.close()
+    vFile.close()
+    rhoFile.close()
+    if params['morphogen']:
+        cFile.close()
+        return (times, mesh, u, v, rho, c)
+    else:
+        return (times, mesh, u, v, rho)
+        
+
+def visualize_1D(params, DIR, interactive=False):
+    if params['morphogen']:
+        (times, mesh, u, v, rho, c) = read_data(params, DIR)
+    else:
+        (times, mesh, u, v, rho) = read_data(params, DIR)
+    u, v = u[:,:,0], v[:,:,0]
+    x = mesh.coordinates()[:,0]
+
+    if params['morphogen']:
+        fig, (axu, axv, axrho, axc) = plt.subplots(1, 4, sharex=True, sharey=True, figsize=(12,3))
+        axc.set_xlabel(r'$x$')
+    else:
+        fig, (axu, axv, axrho) = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(12,4))
+    axu.set_xlabel(r'$x$')
+    axv.set_xlabel(r'$x$')
+    axrho.set_xlabel(r'$x$')
+    axu.set_xlim(np.min(x), np.max(x))
+    axu.set_ylim(np.min(times), np.max(times))
+    axu.set_ylabel(r'$t$')
+    umax, vmax, rhomax = np.max(u), np.max(v), np.max(rho)
+    umin, vmin = -umax, -vmax
+    rhomin = 0
+    ucplot = axu.contourf(x, times, u, cmap='coolwarm', vmin = umin, vmax = umax)
+    vcplot = axv.contourf(x, times, v, cmap='coolwarm', vmin = vmin, vmax = vmax)
+    rhocplot = axrho.contourf(x, times, rho, cmap='viridis', vmin = rhomin, vmax = rhomax)
+    divideru = make_axes_locatable(axu)
+    ucax = divideru.append_axes("right", size="5%", pad=0.05)
+    dividerv = make_axes_locatable(axv)
+    vcax = dividerv.append_axes("right", size="5%", pad=0.05)
+    dividerrho = make_axes_locatable(axrho)
+    rhocax = dividerrho.append_axes("right", size="5%", pad=0.05)
+    plt.colorbar(ucplot, cax=ucax)
+    plt.colorbar(vcplot, cax=vcax)
+    plt.colorbar(rhocplot, cax=rhocax)
+    
+    if params['morphogen']:
+        cmin, cmax = 0, 3*params['average_c']
+        ccplot = axc.contourf(x, times, c, cmap='viridis', vmin = cmin, vmax = cmax)
+        dividerc = make_axes_locatable(axc)
+        ccax = dividerc.append_axes("right", size="5%", pad=0.05)
+        plt.colorbar(ccplot, cax=ccax)
+
+    fig.savefig(os.path.join(DIR, '%s.png' % params['timestamp']), dpi=100)
+    
+    if interactive:
+        if params['morphogen']:
+            fig, (axu, axv, axrho, axc) = plt.subplots(4,1, sharex=True, figsize=(8,8))
+            axc.set_xlabel(r'$x$')
+            axc.set_ylim(np.min(c), np.max(c))
+        else:
+            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])
+        if params['morphogen']:
+            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])
+            if params['morphogen']:
+                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()        
+
+
+def visualize_2D(params, DIR, interactive=False):
+    if params['morphogen']:
+        (times, mesh, u, v, rho, c) = read_data(params, DIR)
+    else:
+        (times, mesh, u, v, rho) = read_data(params, DIR)
+    X = mesh.coordinates()
+    x, y = X[:,0], X[:,1]
+    vx, vy = v[:,:,0], v[:,:,1]
+    rho_min, rho_max = np.min(rho), np.max(rho)
+
+    fig, ax = plt.subplots(1, figsize=(8,8))
+    fig.subplots_adjust(left=0.05, bottom=0.05, top=0.9,
+            right=0.9, wspace=0.0, hspace=0.0)
+    maxspeed = np.sqrt(vx**2+vy**2).max()
+    vx/=maxspeed
+    vy/=maxspeed
+    cax = plt.axes([0.92, 0.05, 0.02, 0.85])
+    levels = np.linspace(rho_min, rho_max, 10)
+    ax.set_aspect(1)
+    ax.axis('off')
+    mesh1 = mesh2triang(mesh)
+    tcf = ax.tricontourf(mesh1, rho[0], levels,
+            cmap='viridis', vmin=rho_min, vmax=rho_max, antialiased=True)
+    cbar = fig.colorbar(tcf, cax=cax)
+    cbar.set_label(r'$c/c_o$', labelpad=-10)
+    lev_labels = [levels[0], levels[-1]]
+    cbar.set_ticks(lev_labels)
+    cbar.ax.set_yticklabels(['{:3.1f}'.format(x) for x in lev_labels])
+    qui = ax.quiver(x, y, vx[0], vy[0],
+            edgecolors='k', color='k', pivot='mid',
+            scale=32, linewidth=0.1, headwidth=3)
+    def update2(val):
+        ti = (abs(times-val)).argmin()
+        ax.clear()
+        ax.axis('off')
+        ax.tricontourf(mesh1, rho[ti], levels,
+            cmap='viridis', vmin=rho_min, vmax=rho_max, antialiased=True)
+        ax.quiver(x, y, vx[ti], vy[ti], 
+            edgecolors='k', color='k', pivot='mid',
+            scale=32, linewidth=0.1, headwidth=3)
+        plt.draw()
+    sax = plt.axes([0.05, 0.92, 0.85, 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(update2)
+
+    if interactive:
+        plt.show()
+    else:
+        print('Saving movie...')
+        FPS = 10
+        movFile = os.path.join(DIR, '%s.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))):
+            idx = (abs(times-times[tt])).argmin()
+            update2(idx)
+            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()
+
+
+def visualize(params, DIR, interactive=False):
+    if params['dimension']==1:
+        visualize_1D(params, DIR, interactive)
+    elif params['dimension']==2:
+        visualize_2D(params, DIR, interactive)
+        
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser()
+    parser.add_argument('-j','--jsonfile', help='json file', default='fixed_boundaries.json')
+    args = parser.parse_args()
+
+    with open(args.jsonfile) as jsonFile:
+        params = json.load(jsonFile)
+    DIR = os.path.dirname(args.jsonfile)
+    
+    visualize(params, DIR, interactive=True)
+
+
+
+
+
+