fixed bug in calling visualize, DIR

fixed_boundaries.json

 {
-    "morphogen" : true,
+    "morphogen" : false,
 
-    "resolution" : 64,
+    "resolution" : 32,
     "system_size_by_2PI" : 1 ,
 
     "elasticity" : 0.0,
     "viscosity" : 1.0,
     "friction" : 1.0,
-    "lamda": 1.75,
+    "lamda": 5.5,
 
     "diffusion_rho" : 1.0,
     "turnover_rho" : 0.0,
     "average_rho" : 1.0,
     "saturation_rho" : 1.0,
 
-    "diffusion_c" : 1.0,
+    "diffusion_c" : 0.0,
     "turnover_c" : 0.0,
-    "average_c" : 1.0,
+    "average_c" : 0.0,
 
     "noise_level": 0.1,
 
-    "timestep" : 0.05,
-    "maxtime" : 10.0
+    "timestep" : 0.01,
+    "savetime": 0.5,
+    "maxtime" : 200.0
 
 }

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
@@ -156,13 +157,13 @@ class FixedBoundaries(object):
         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):
+    def solve(self, extend=False, DIR=''):
         # for saving data
-        self.uFile   = df.XDMFFile('%s_displacement.xdmf' % self.timestamp)
-        self.vFile   = df.XDMFFile('%s_velocity.xdmf' % self.timestamp)
-        self.rhoFile = df.XDMFFile('%s_density.xdmf' % self.timestamp)
+        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('%s_concentration.xdmf' % self.timestamp)
+            self.cFile = df.XDMFFile(os.path.join(DIR, '%s_concentration.xdmf' % self.timestamp))
         
         if extend:
             SFS = df.FunctionSpace(self.mesh, 'P', 1)
@@ -180,13 +181,19 @@ class FixedBoundaries(object):
         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 i in progressbar.progressbar(range(1, int(self.maxtime/self.timestep))):
+        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()

run_fixed_boundaries.py

@@ -4,6 +4,7 @@ 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')
@@ -33,7 +34,7 @@ else:
 
 # solve
 fb = FixedBoundaries(params)
-fb.solve(extend=extend)
+fb.solve(extend=extend, DIR=DIR)
 
 if extend:
     params['maxtime'] = oldMaxTime + params['maxtime']
@@ -41,3 +42,5 @@ if extend:
 with open(os.path.join(DIR, 
             params['timestamp'] + '_fixed_boundaries.json'), "w") as fp:
     json.dump(params, fp, indent=4)
+
+visualize(params, DIR, interactive=False)
\ No newline at end of file

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" : True,
+    "morphogen" : False,
 
     "resolution" : 32,
     "system_size_by_2PI" : 1 ,
@@ -17,37 +21,39 @@ parameters = {
     "elasticity" : 0.0,
     "viscosity" : 1.0,
     "friction" : 1.0,
-    "lamda": 1.75,
+    "lamda": 5.5,
 
     "diffusion_rho" : 1.0,
     "turnover_rho" : 0.0,
     "average_rho" : 1.0,
     "saturation_rho" : 1.0,
 
-    "diffusion_c" : 1.0,
+    "diffusion_c" : 0.0,
     "turnover_c" : 0.0,
-    "average_c" : 1.0,
+    "average_c" : 0.0,
 
     "noise_level": 0.1,
 
-    "timestep" : 0.05,
-    "maxtime" : 10.0
+    "timestep" : 0.01,
+    "savetime": 0.5,
+    "maxtime" : 200.0
 
 }
 
 # parameters to scan
-lamda_list = np.arange(0, 2.2, 0.2)
-elasticity_list = np.arange(0, 1.1, 0.1)
+lamda_list = np.arange(1.0, 61.0, 10.0)
+turnover_rho_list = np.arange(0, 1.0, 0.01)
+elasticity_list = np.arange(0, 1.0, 0.01)
 
-plist = list(itertools.product(lamda_list, elasticity_list))
+plist = list(itertools.product(lamda_list, turnover_rho_list, elasticity_list))
 
 # scan parameters
 for i in range(len(plist)):
-    lamda, elasticity = plist[i]
+    lamda, turnover_rho, elasticity = plist[i]
 
     print('-'*50)
     print('%d/%d' % (i, len(plist)))
-    print('lamda = %4.3f, elasticity = %4.3f' % (lamda, elasticity))
+    print('lamda = %4.3f, turnover_rho = %4.3f, elasticity = %4.3f' % (lamda, turnover_rho, elasticity))
     print('-'*50)
 
     # timestamp as job-ID
@@ -55,15 +61,17 @@ for i in range(len(plist)):
 
     params = parameters.copy()
     params['timestamp'] = timestamp
-    params['lamda'] = lamda
+    params['turnover_rho'] = turnover_rho
     params['elasticity'] = elasticity
 
     try:
         fb = FixedBoundaries(params)
         fb.solve(extend=False)
-
-        with open(timestamp+'_polarity.json', "w") as fp:
+        DIR = args.outputdir
+        with open(os.path.join(DIR, params['timestamp']+'_fixed_boundaries.json'), "w") as fp:
             json.dump(params, fp, indent=4)
+            
+        visualize(params, DIR, interactive=False)
     except:
             print('*** FAILED ****')
             os.system("rm -f %s*" % timestamp)

test

-yyy
\ No newline at end of file

viz_fixed_boundaries.py

@@ -3,33 +3,30 @@ import json
 import matplotlib.pyplot as plt
 from matplotlib.widgets import Slider
 import argparse
+import progressbar
+import os
 import dolfin as df
 
-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()
+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']
     
-with open(args.jsonfile) as jsonFile:
-    params = json.load(jsonFile)
-
-mesh = df.IntervalMesh(params['resolution'], 0.0, 2.0*np.pi*params['system_size_by_2PI'])
-x = mesh.coordinates()
-SFS = df.FunctionSpace(mesh, 'P', 1)
-ui, vi, rhoi, ci = df.Function(SFS), df.Function(SFS), df.Function(SFS), df.Function(SFS)
+    mesh = df.IntervalMesh(params['resolution'], 0.0, 2.0*np.pi*params['system_size_by_2PI'])
+    x = mesh.coordinates()[:,0]
+    SFS = df.FunctionSpace(mesh, 'P', 1)
+    ui, vi, rhoi = df.Function(SFS), df.Function(SFS), df.Function(SFS)
 
-savesteps = int(params['maxtime']/params['timestep'])
-times =  np.arange(savesteps+1) * params['timestep']
-
-u = np.zeros((len(times), len(x)))
-v, rho = np.zeros_like(u), np.zeros_like(u)
-uFile   = df.XDMFFile('%s_displacement.xdmf' % params['timestamp'])
-vFile   = df.XDMFFile('%s_velocity.xdmf' % params['timestamp'])
-rhoFile = df.XDMFFile('%s_density.xdmf' % params['timestamp'])
-if params['morphogen']:
+    u = np.zeros((len(times), len(x)))
+    v, rho = np.zeros_like(u), np.zeros_like(u)
+    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(u)
-    cFile = df.XDMFFile('%s_concentration.xdmf' % params['timestamp'])
-for steps in range(len(times)):
+        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)
@@ -39,33 +36,65 @@ for steps in range(len(times)):
         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']:
+    uFile.close()
+    vFile.close()
+    rhoFile.close()
+    if params['morphogen']:
         cFile.close()
+        return (times, x, u, v, rho, c)
+    else:
+        return (times, x, u, v, rho)
+        
+
+def visualize(params, DIR, interactive=False):
+    if params['morphogen']:
+        (times, x, u, v, rho, c) = read_data(params, DIR)
+    else:
+        (times, x, u, v, rho) = read_data(params, DIR)
+
+    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$')
+    
+    axu.contourf(x, times, u, cmap='coolwarm')
+    axv.contourf(x, times, v, cmap='coolwarm')
+    axrho.contourf(x, times, rho, cmap='viridis')
+    
+    if params['morphogen']:
+        axc.contourf(x, times, c, cmap='viridis')
         
-if params['morphogen']:
+    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:
+        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']:
+        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):
+        def update(value):
             ti = np.abs(times-value).argmin()
             uplot.set_ydata(u[ti])
             vplot.set_ydata(v[ti])
@@ -74,10 +103,27 @@ def update(value):
                 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),
+        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()
+        slider.drawon = False
+        slider.on_changed(update)
+        plt.show()        
+
+
+
+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)
+
+
+
+