Skip to content

Jigyasa Watwani / growth-pattern-control

Go to a project
Commit 4dd8bff7 by Vijay Kumar Krishnamurthy
Options

Fully working code. Fixed bugs

parent c9aa2b3e
Inline Side-by-side
Showing with 121 additions and 60 deletions
fixed_boundaries.json
{ {
"morphogen" : true, "morphogen" : false,
"resolution" : 64, "resolution" : 32,
"system_size_by_2PI" : 1 , "system_size_by_2PI" : 1 ,
"elasticity" : 0.0, "elasticity" : 0.0,
"viscosity" : 1.0, "viscosity" : 1.0,
"friction" : 1.0, "friction" : 1.0,
"lamda": 1.75, "lamda": 5.5,
"diffusion_rho" : 1.0, "diffusion_rho" : 1.0,
"turnover_rho" : 0.0, "turnover_rho" : 0.0,
"average_rho" : 1.0, "average_rho" : 1.0,
"saturation_rho" : 1.0, "saturation_rho" : 1.0,
"diffusion_c" : 1.0, "diffusion_c" : 0.0,
"turnover_c" : 0.0, "turnover_c" : 0.0,
"average_c" : 1.0, "average_c" : 0.0,
"noise_level": 0.1, "noise_level": 0.1,
"timestep" : 0.05, "timestep" : 0.01,
"maxtime" : 10.0 "savetime": 0.5,
"maxtime" : 200.0
} }
fixed_boundaries.py
import dolfin as df import dolfin as df
import numpy as np import numpy as np
import progressbar import progressbar
import os
df.set_log_level(df.LogLevel.ERROR) df.set_log_level(df.LogLevel.ERROR)
df.parameters['form_compiler']['optimize'] = True df.parameters['form_compiler']['optimize'] = True
...@@ -156,13 +157,13 @@ class FixedBoundaries(object): ...@@ -156,13 +157,13 @@ class FixedBoundaries(object):
self.rhoFile.write_checkpoint(rho, 'density', self.time, append=True) self.rhoFile.write_checkpoint(rho, 'density', self.time, append=True)
self.cFile.write_checkpoint(c, 'concentration', 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 # for saving data
self.uFile = df.XDMFFile('%s_displacement.xdmf' % self.timestamp) self.uFile = df.XDMFFile(os.path.join(DIR, '%s_displacement.xdmf' % self.timestamp))
self.vFile = df.XDMFFile('%s_velocity.xdmf' % self.timestamp) self.vFile = df.XDMFFile(os.path.join(DIR, '%s_velocity.xdmf' % self.timestamp))
self.rhoFile = df.XDMFFile('%s_density.xdmf' % self.timestamp) self.rhoFile = df.XDMFFile(os.path.join(DIR, '%s_density.xdmf' % self.timestamp))
if self.morphogen: 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: if extend:
SFS = df.FunctionSpace(self.mesh, 'P', 1) SFS = df.FunctionSpace(self.mesh, 'P', 1)
...@@ -180,13 +181,19 @@ class FixedBoundaries(object): ...@@ -180,13 +181,19 @@ class FixedBoundaries(object):
self.setup_initial_conditions(ui, rhoi, ci) self.setup_initial_conditions(ui, rhoi, ci)
self.setup_weak_forms() self.setup_weak_forms()
# time-variables
self.time = 0.0 self.time = 0.0
savesteps = int(self.savetime/self.timestep)
maxsteps = int(self.maxtime/self.timestep)
if not extend:
self.savedata() 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) df.solve(self.form == 0, self.f, self.bc)
self.time += self.timestep self.time += self.timestep
df.assign(self.f1, self.f0) df.assign(self.f1, self.f0)
df.assign(self.f0, self.f) df.assign(self.f0, self.f)
if steps % savesteps == 0:
self.savedata() self.savedata()
self.uFile.close() self.uFile.close()
......
run_fixed_boundaries.py
...@@ -4,6 +4,7 @@ import os ...@@ -4,6 +4,7 @@ import os
import argparse import argparse
from fixed_boundaries import FixedBoundaries from fixed_boundaries import FixedBoundaries
from viz_fixed_boundaries import visualize
parser = argparse.ArgumentParser() parser = argparse.ArgumentParser()
parser.add_argument('-j','--jsonfile', help='json file', default='fixed_boundaries.json') parser.add_argument('-j','--jsonfile', help='json file', default='fixed_boundaries.json')
...@@ -33,7 +34,7 @@ else: ...@@ -33,7 +34,7 @@ else:
# solve # solve
fb = FixedBoundaries(params) fb = FixedBoundaries(params)
fb.solve(extend=extend) fb.solve(extend=extend, DIR=DIR)
if extend: if extend:
params['maxtime'] = oldMaxTime + params['maxtime'] params['maxtime'] = oldMaxTime + params['maxtime']
...@@ -41,3 +42,5 @@ if extend: ...@@ -41,3 +42,5 @@ if extend:
with open(os.path.join(DIR, with open(os.path.join(DIR,
params['timestamp'] + '_fixed_boundaries.json'), "w") as fp: params['timestamp'] + '_fixed_boundaries.json'), "w") as fp:
json.dump(params, fp, indent=4) json.dump(params, fp, indent=4)
visualize(params, DIR, interactive=False)
\ No newline at end of file
scan_fixed_boundaries.py
...@@ -3,6 +3,7 @@ import dolfin as df ...@@ -3,6 +3,7 @@ import dolfin as df
import numpy as np import numpy as np
from fixed_boundaries import FixedBoundaries from fixed_boundaries import FixedBoundaries
from viz_fixed_boundaries import visualize
df.set_log_level(df.LogLevel.ERROR) df.set_log_level(df.LogLevel.ERROR)
df.parameters['form_compiler']['optimize'] = True df.parameters['form_compiler']['optimize'] = True
...@@ -64,6 +65,8 @@ for i in range(len(plist)): ...@@ -64,6 +65,8 @@ for i in range(len(plist)):
with open(timestamp+'_polarity.json', "w") as fp: with open(timestamp+'_polarity.json', "w") as fp:
json.dump(params, fp, indent=4) json.dump(params, fp, indent=4)
visualize(params, interactive=False)
except: except:
print('*** FAILED ****') print('*** FAILED ****')
os.system("rm -f %s*" % timestamp) os.system("rm -f %s*" % timestamp)
......
viz_fixed_boundaries.py
...@@ -3,33 +3,30 @@ import json ...@@ -3,33 +3,30 @@ import json
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from matplotlib.widgets import Slider from matplotlib.widgets import Slider
import argparse import argparse
import progressbar
import os
import dolfin as df import dolfin as df
parser = argparse.ArgumentParser() def read_data(params, DIR):
parser.add_argument('-j','--jsonfile', help='json file', default='fixed_boundaries.json') print('Reading data from %s/%s...' % (DIR, params['timestamp']))
parser.add_argument('-o','--outputdir', help='output directory', type=str, default='') savesteps = int(params['maxtime']/params['savetime'])
args = parser.parse_args() times = np.arange(savesteps+1) * params['savetime']
with open(args.jsonfile) as jsonFile: mesh = df.IntervalMesh(params['resolution'], 0.0, 2.0*np.pi*params['system_size_by_2PI'])
params = json.load(jsonFile) x = mesh.coordinates()[:,0]
SFS = df.FunctionSpace(mesh, 'P', 1)
mesh = df.IntervalMesh(params['resolution'], 0.0, 2.0*np.pi*params['system_size_by_2PI']) ui, vi, rhoi = df.Function(SFS), df.Function(SFS), df.Function(SFS)
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)
savesteps = int(params['maxtime']/params['timestep']) u = np.zeros((len(times), len(x)))
times = np.arange(savesteps+1) * params['timestep'] v, rho = np.zeros_like(u), np.zeros_like(u)
uFile = df.XDMFFile(os.path.join(DIR, '%s_displacement.xdmf' % params['timestamp']))
u = np.zeros((len(times), len(x))) vFile = df.XDMFFile(os.path.join(DIR, '%s_velocity.xdmf' % params['timestamp']))
v, rho = np.zeros_like(u), np.zeros_like(u) rhoFile = df.XDMFFile(os.path.join(DIR, '%s_density.xdmf' % params['timestamp']))
uFile = df.XDMFFile('%s_displacement.xdmf' % params['timestamp']) if params['morphogen']:
vFile = df.XDMFFile('%s_velocity.xdmf' % params['timestamp']) ci = df.Function(SFS)
rhoFile = df.XDMFFile('%s_density.xdmf' % params['timestamp'])
if params['morphogen']:
c = np.zeros_like(u) c = np.zeros_like(u)
cFile = df.XDMFFile('%s_concentration.xdmf' % params['timestamp']) cFile = df.XDMFFile(os.path.join(DIR, '%s_concentration.xdmf' % params['timestamp']))
for steps in range(len(times)): for steps in progressbar.progressbar(range(savesteps+1)):
uFile.read_checkpoint(ui, 'displacement', steps) uFile.read_checkpoint(ui, 'displacement', steps)
vFile.read_checkpoint(vi, 'velocity', steps) vFile.read_checkpoint(vi, 'velocity', steps)
rhoFile.read_checkpoint(rhoi, 'density', steps) rhoFile.read_checkpoint(rhoi, 'density', steps)
...@@ -39,33 +36,65 @@ for steps in range(len(times)): ...@@ -39,33 +36,65 @@ for steps in range(len(times)):
if params['morphogen']: if params['morphogen']:
cFile.read_checkpoint(ci, 'concentration', steps) cFile.read_checkpoint(ci, 'concentration', steps)
c[steps] = ci.compute_vertex_values(mesh) c[steps] = ci.compute_vertex_values(mesh)
uFile.close() uFile.close()
vFile.close() vFile.close()
rhoFile.close() rhoFile.close()
if params['morphogen']: if params['morphogen']:
cFile.close() 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)) fig, (axu, axv, axrho, axc) = plt.subplots(4,1, sharex=True, figsize=(8,8))
axc.set_xlabel(r'$x$') axc.set_xlabel(r'$x$')
axc.set_ylim(np.min(c), np.max(c)) axc.set_ylim(np.min(c), np.max(c))
else: else:
fig, (axu, axv, axrho) = plt.subplots(3,1, sharex=True, figsize=(8,6)) fig, (axu, axv, axrho) = plt.subplots(3,1, sharex=True, figsize=(8,6))
axrho.set_xlabel(r'$x$') axrho.set_xlabel(r'$x$')
axu.set_ylabel(r'$u(x,t)$') axu.set_ylabel(r'$u(x,t)$')
axv.set_ylabel(r'$v(x,t)$') axv.set_ylabel(r'$v(x,t)$')
axrho.set_ylabel(r'$\rho(x,t)$') axrho.set_ylabel(r'$\rho(x,t)$')
axu.set_ylim(np.min(u), np.max(u)) axu.set_ylim(np.min(u), np.max(u))
axv.set_ylim(np.min(v), np.max(v)) axv.set_ylim(np.min(v), np.max(v))
axrho.set_ylim(np.min(rho), np.max(rho)) axrho.set_ylim(np.min(rho), np.max(rho))
uplot, = axu.plot(x, u[0]) uplot, = axu.plot(x, u[0])
vplot, = axv.plot(x, v[0]) vplot, = axv.plot(x, v[0])
rhoplot, = axrho.plot(x, rho[0]) rhoplot, = axrho.plot(x, rho[0])
if params['morphogen']: if params['morphogen']:
cplot, = axc.plot(x, c[0]) cplot, = axc.plot(x, c[0])
def update(value): def update(value):
ti = np.abs(times-value).argmin() ti = np.abs(times-value).argmin()
uplot.set_ydata(u[ti]) uplot.set_ydata(u[ti])
vplot.set_ydata(v[ti]) vplot.set_ydata(v[ti])
...@@ -74,10 +103,27 @@ def update(value): ...@@ -74,10 +103,27 @@ def update(value):
cplot.set_ydata(c[ti]) cplot.set_ydata(c[ti])
plt.draw() plt.draw()
sax = plt.axes([0.1, 0.92, 0.7, 0.02]) sax = plt.axes([0.1, 0.92, 0.7, 0.02])
slider = Slider(sax, r'$t/\tau$', min(times), max(times), slider = Slider(sax, r'$t/\tau$', min(times), max(times),
valinit=min(times), valfmt='%3.1f', valinit=min(times), valfmt='%3.1f',
fc='#999999') fc='#999999')
slider.drawon = False slider.drawon = False
slider.on_changed(update) slider.on_changed(update)
plt.show() 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)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or sign in to comment