extend

fixed_boundaries.json

@@ -21,6 +21,6 @@
     "noise_level": 0.1,
     "noise_level": 0.1,
 
 
     "timestep" : 0.05,
     "timestep" : 0.05,
-    "maxtime" : 200.0
+    "maxtime" : 10.0
 
 
 }
 }

fixed_boundaries.py

@@ -60,8 +60,18 @@ class FixedBoundaries(object):
     def diffusion_reaction_c(self, c, tc):
     def diffusion_reaction_c(self, c, tc):
         return (self.diffusion_c * df.inner(c.dx(0), tc.dx(0)) + self.reaction_c(c, tc))
         return (self.diffusion_c * df.inner(c.dx(0), tc.dx(0)) + self.reaction_c(c, tc))
 
 
-    def setup_initial_conditions(self, u0=None, rho0=None, c0=None):
-        if u0==None and rho0==None and c0==None:
+    def setup_initial_conditions(self, ui=None, rhoi=None, ci=None):
+        if (ui is not None) and (rhoi is not None):
+            u0   = df.Function(self.function_space.sub(0).collapse())
+            u0.vector()[:] = ui
+            rho0 = df.Function(self.function_space.sub(2).collapse())
+            rho0.vector()[:] = rhoi
+            if self.morphogen:
+                c0 = df.Function(self.function_space.sub(3).collapse())
+                c0.vector()[:] = ci
+            else:
+                c0 = df.Constant(1.0)
+        else:
             u0   = df.interpolate(df.Constant(0), self.function_space.sub(0).collapse())
             u0   = df.interpolate(df.Constant(0), self.function_space.sub(0).collapse())
             rho0 = df.interpolate(df.Constant(self.average_rho), self.function_space.sub(2).collapse())
             rho0 = df.interpolate(df.Constant(self.average_rho), self.function_space.sub(2).collapse())
             # add noise
             # add noise
@@ -76,13 +86,7 @@ class FixedBoundaries(object):
                 c0.vector()[:] = c0.vector()[:] + noise_c
                 c0.vector()[:] = c0.vector()[:] + noise_c
             else:
             else:
                 c0 = df.Constant(1.0)
                 c0 = df.Constant(1.0)
-        else:
-            u0   = df.interpolate(u0, self.function_space.sub(0).collapse())
-            rho0 = df.interpolate(rho0, self.function_space.sub(2).collapse())
-            if self.morphogen:
-                c0 = df.interpolate(c0, self.function_space.sub(3).collapse())
-            else:
-                c0 = df.Constant(1.0)
+
                 
                 
         VFS = self.function_space.sub(1).collapse()
         VFS = self.function_space.sub(1).collapse()
         v0  = df.Function(VFS)
         v0  = df.Function(VFS)
@@ -145,13 +149,7 @@ class FixedBoundaries(object):
         Z = f.split(deepcopy=True)
         Z = f.split(deepcopy=True)
         return np.array([z.compute_vertex_values(self.mesh) for z in Z]).T
         return np.array([z.compute_vertex_values(self.mesh) for z in Z]).T
 
 
-    def solve(self, extend=None, DIR=''):
-        if extend:
-            # read the last time point of earlier simulation
-            # u0, rho0, c0
-            pass
-        else:
-            u0, rho0, c0 = None, None, None
+    def solve(self, u0=None, rho0=None, c0=None):
         self.setup_initial_conditions(u0, rho0, c0)
         self.setup_initial_conditions(u0, rho0, c0)
         self.setup_weak_forms()
         self.setup_weak_forms()
 
 
@@ -165,69 +163,5 @@ class FixedBoundaries(object):
             data.append(self.get_data(self.f))
             data.append(self.get_data(self.f))
             self.f0.assign(self.f)
             self.f0.assign(self.f)
             self.f1.assign(self.f0)
             self.f1.assign(self.f0)
-            
-        # parse data
-        data = np.array(data)
-        if self.morphogen:
-            # u, v, rho, c
-            return np.split(data, 4, axis=2)
-        else:
-            # u, v, rho
-            return np.split(data, 3, axis=2)
-
 
 
-if __name__ == '__main__':
-    import json
+        return np.array(data)
\ No newline at end of file
-    import matplotlib.pyplot as plt
-    from matplotlib.widgets import Slider
-    
-    with open('fixed_boundaries.json') as jsonFile:
-        params = json.load(jsonFile)
-
-    fb = FixedBoundaries(params)
-    Z = fb.solve()
-    
-    if params['morphogen']:
-        u, v, rho, c = Z
-    else:
-        u, v, rho = Z
-
-
-    x = fb.mesh.coordinates()
-    times = np.arange(0, fb.maxtime, fb.timestep)
-    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()

run_fixed_boundaries.py

+import datetime
+import json
+import os
+import argparse
+
+import numpy as np
+
+from fixed_boundaries import FixedBoundaries
+
+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
+    dataFile = params['timestamp'] + '_fixed_boundaries.npy'
+    assert os.path.isfile(dataFile), '%s not found' % dataFile
+    old_data = np.load(dataFile)
+    print(old_data.shape)
+    u0, rho0 = old_data[-1, :, 0], old_data[-1, :, 2]
+    if params['morphogen']:
+        c0 = old_data[-1, :, 3]
+
+# solve
+fb = FixedBoundaries(params)
+data = fb.solve(u0=u0, rho0=rho0, c0=c0)
+
+if extend:
+    params['maxtime'] = oldMaxTime + params['maxtime']
+    data = np.concatenate((old_data, data[1:]), axis=0)
+
+# save data
+np.save('%s_fixed_boundaries.npy' % params['timestamp'], data)
+        
+with open(os.path.join(DIR, 
+            params['timestamp'] + '_fixed_boundaries.json'), "w") as fp:
+    json.dump(params, fp, indent=4)
\ No newline at end of file

viz_fixed_boundaries.py

+import numpy as np
+import json
+import matplotlib.pyplot as plt
+from matplotlib.widgets import Slider
+
+import argparse
+
+parser = argparse.ArgumentParser()
+parser.add_argument('-j','--jsonfile', help='json file', default='fixed_boundaries.json')
+parser.add_argument('-o','--outputdir', help='output directory', type=str, default='')
+args = parser.parse_args()
+
+with open(args.jsonfile) as jsonFile:
+    params = json.load(jsonFile)
+
+dataFile = params['timestamp'] + '_fixed_boundaries.npy'
+data = np.load(dataFile)
+
+u, v, rho = data[:,:,0], data[:,:,1], data[:,:,2]
+if params['morphogen']:
+    c = data[:,:,3]
+
+
+x = np.linspace(0, params['system_size_by_2PI']*2*np.pi, params['resolution']+1)
+times = np.arange(0, params['maxtime'], params['timestep'])
+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()