test

fixed_boundaries.py

@@ -18,9 +18,11 @@ class FixedBoundaries(object):
             # u, v, rho, c
             mixed_element = df.MixedElement([scalar_element, scalar_element, 
                                              scalar_element, scalar_element]) 
+            self.savedata = self.save_data_uvrhoc
         else:
             # u, v, rho
             mixed_element = df.MixedElement([scalar_element, scalar_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)
@@ -62,13 +64,10 @@ class FixedBoundaries(object):
 
     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
+            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.Function(self.function_space.sub(3).collapse())
-                c0.vector()[:] = ci
+                c0 = df.interpolate(ci, self.function_space.sub(3).collapse())
             else:
                 c0 = df.Constant(1.0)            
         else:
@@ -87,7 +86,6 @@ class FixedBoundaries(object):
             else:
                 c0 = df.Constant(1.0)
         
-                
         VFS = self.function_space.sub(1).collapse()
         v0  = df.Function(VFS)
         tv = df.TestFunction(VFS)
@@ -145,23 +143,54 @@ class FixedBoundaries(object):
         else:
             self.form = (uform + vform + rhoform) * df.dx
             
-    def get_data(self, f):
-        Z = f.split(deepcopy=True)
-        return np.array([z.compute_vertex_values(self.mesh) for z in Z]).T
+    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):
+        # 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)
+        if self.morphogen:
+            self.cFile = df.XDMFFile('%s_concentration.xdmf' % self.timestamp)
+        
+        if extend:
+            SFS = df.FunctionSpace(self.mesh, 'P', 1)
+            ui, vi, rhoi, ci = df.Function(SFS), df.Function(SFS), 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     
             
-    def solve(self, u0=None, rho0=None, c0=None):
-        self.setup_initial_conditions(u0, rho0, c0)
+        self.setup_initial_conditions(ui, rhoi, ci)
         self.setup_weak_forms()
         
-        # store
-        data = []
-        data.append(self.get_data(self.f0))
-        
-        times = np.arange(0.0, self.maxtime, self.timestep)
-        for i in progressbar.progressbar(range(0, len(times))):
+        self.time = 0.0
+        self.savedata()
+        for i in progressbar.progressbar(range(1, int(self.maxtime/self.timestep))):
             df.solve(self.form == 0, self.f, self.bc)
-            data.append(self.get_data(self.f))
-            self.f0.assign(self.f)
-            self.f1.assign(self.f0)
-
-        return np.array(data)
\ No newline at end of file
+            self.time += self.timestep
+            df.assign(self.f1, self.f0)
+            df.assign(self.f0, self.f)  
+            self.savedata()
+
+        self.uFile.close()
+        self.vFile.close()
+        self.rhoFile.close()
+        if self.morphogen:
+            self.cFile.close()
\ No newline at end of file

run_fixed_boundaries.py

@@ -3,8 +3,6 @@ import json
 import os
 import argparse
 
-import numpy as np
-
 from fixed_boundaries import FixedBoundaries
 
 parser = argparse.ArgumentParser()
@@ -32,24 +30,13 @@ else:
     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)
+fb.solve(extend=extend)
 
 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:

scan_fixed_boundaries.py

+import datetime, time, json, os, itertools
+import dolfin as df
+import numpy as np
+
+from fixed_boundaries import FixedBoundaries
+
+df.set_log_level(df.LogLevel.ERROR)
+df.parameters['form_compiler']['optimize'] = True
+
+# base parameters
+parameters = {
+    "morphogen" : True,
+
+    "resolution" : 32,
+    "system_size_by_2PI" : 1 ,
+
+    "elasticity" : 0.0,
+    "viscosity" : 1.0,
+    "friction" : 1.0,
+    "lamda": 1.75,
+
+    "diffusion_rho" : 1.0,
+    "turnover_rho" : 0.0,
+    "average_rho" : 1.0,
+    "saturation_rho" : 1.0,
+
+    "diffusion_c" : 1.0,
+    "turnover_c" : 0.0,
+    "average_c" : 1.0,
+
+    "noise_level": 0.1,
+
+    "timestep" : 0.05,
+    "maxtime" : 10.0
+
+}
+
+# parameters to scan
+lamda_list = np.arange(0, 2.2, 0.2)
+elasticity_list = np.arange(0, 1.1, 0.1)
+
+plist = list(itertools.product(lamda_list, elasticity_list))
+
+# scan parameters
+for i in range(len(plist)):
+    lamda, elasticity = plist[i]
+
+    print('-'*50)
+    print('%d/%d' % (i, len(plist)))
+    print('lamda = %4.3f, elasticity = %4.3f' % (lamda, 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['lamda'] = lamda
+    params['elasticity'] = elasticity
+
+    try:
+        fb = FixedBoundaries(params)
+        fb.solve(extend=False)
+
+        with open(timestamp+'_polarity.json', "w") as fp:
+            json.dump(params, fp, indent=4)
+    except:
+        print('*** FAILED ****')
+        os.system("rm -f %s*" % timestamp)
+
+    time.sleep(2)

test

+yyy
\ No newline at end of file

viz_fixed_boundaries.py

@@ -2,8 +2,8 @@ import numpy as np
 import json
 import matplotlib.pyplot as plt
 from matplotlib.widgets import Slider
-
 import argparse
+import dolfin as df
 
 parser = argparse.ArgumentParser()
 parser.add_argument('-j','--jsonfile', help='json file', default='fixed_boundaries.json')
@@ -13,16 +13,38 @@ args = parser.parse_args()
 with open(args.jsonfile) as jsonFile:
     params = json.load(jsonFile)
 
-dataFile = params['timestamp'] + '_fixed_boundaries.npy'
-data = np.load(dataFile)
+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)
             
-u, v, rho = data[:,:,0], data[:,:,1], data[:,:,2]
-if params['morphogen']:
-    c = data[:,:,3]
+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']:
+    c = np.zeros_like(u)
+    cFile = df.XDMFFile('%s_concentration.xdmf' % params['timestamp'])
+for steps in range(len(times)):
+    uFile.read_checkpoint(ui, 'displacement', steps)
+    vFile.read_checkpoint(vi, 'velocity', steps)
+    rhoFile.read_checkpoint(rhoi, 'density', steps)
+    u[steps] = ui.compute_vertex_values(mesh)
+    v[steps] = vi.compute_vertex_values(mesh)
+    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()
 
-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$')