Skip to content

Jigyasa Watwani / growth-pattern-control

Go to a project
Commit 9e0de96a by Jigyasa Watwani
Options

separated visualisation, code to scan parameters and plot phase diagram

parent be194ebf
Inline Side-by-side
Showing with 275 additions and 85 deletions
basic_model_only_density/phase_diagram_tissue.py 0 → 100644
#!/usr/bin/env python3
import dolfin as df
import numpy as np
import glob
import json
import h5py
import progressbar
import os
import argparse
import matplotlib.pyplot as plt
from matplotlib.backend_bases import MouseButton
from scipy import spatial
def phaseDiagram(DIR):
J = glob.glob(DIR+'/*.json') # list of files with this path
# J.remove(DIR+'/parameters.json')
parameters = dict()
peclet = []
retardationtime = []
timeStamper = []
colList = []
# check data files
for j in J:
timestamp = os.path.basename(j).split('_')[0]
dataFile = ['%s_density.xdmf' % timestamp,
'%s_density.h5' % timestamp,
'%s_velocity.xdmf' % timestamp,
'%s_velocity.h5' % timestamp,
'%s_displacement.xdmf' % timestamp,
'%s_displacement.h5' % timestamp]
for d in dataFile:
assert os.path.isfile(os.path.join(DIR,d)), '%s not found' % d
with open(j) as jFile:
try:
p = json.load(jFile)
# del p['timestamp']
try:
length = np.load('%s_length.npy' % p['timestamp'])
if (length[-1]-length[-3]<1e-7):
print('controlled_growth')
col='#e9c46a'
else:
print('uncontrolled_growth')
col='#2a9d8f'
colList.append(col)
parameters[timestamp] = p #?
Pe = (p['lamda'] /
(p['viscosity']*p['turnover_rho']))
rt = (p['elasticity']
/ (p['viscosity']*p['turnover_rho']))
peclet.append(Pe)
retardationtime.append(rt)
timeStamper.append(timestamp)
except:
print('Unable to read %s_length.npy' % timestamp)
except:
print('Error opening %s' % j)
fig, ax = plt.subplots(1)
for i in range(len(retardationtime)):
ax.scatter(retardationtime[i], peclet[i], s=100, fc=colList[i])
ax.set_xlabel(r'$K/k_\rho \eta$')
ax.set_ylabel(r'$\lambda /k_\rho \eta$')
plt.show()
if __name__=="__main__":
import argparse
parser = argparse.ArgumentParser()
parser.add_argument('-d','--directory',
help='directory with json files', required=True)
args = parser.parse_args()
assert os.path.isdir(args.directory), \
'%s directory not found' % args.directory
phaseDiagram(args.directory)
basic_model_only_density/scan_tissue.py 0 → 100644
from __future__ import print_function
import json
import itertools
import os
import datetime
import time
import numpy as np
from tissue import Tissue
assert os.path.isfile("parameters.json"), 'parameters.json file not found'
with open("parameters.json") as jsonFile:
parameters = json.load(jsonFile)
#=========================================
# parameters to scan
lamda = np.arange(-5, -0.0, 0.5)
elasticity = np.arange(0.2, 5.2, 0.5)
plist = list(itertools.product(lamda, elasticity))
#=========================================
# scan
for p in plist:
lamda, elasticity = p
print('\n------------------\n')
print('lamda = ',lamda, ', elasticity = ', elasticity)
timestamp = datetime.datetime.now().strftime("%d%m%y-%H%M%S")
params = parameters.copy()
params['timestamp'] = timestamp
params['lamda'] = lamda
params['elasticity'] = elasticity
t = Tissue(params)
t.solve()
parametersFile = timestamp + '__parameters.json'
with open(parametersFile, "w") as fp:
json.dump(params, fp, indent=4, sort_keys=True, default=str)
from viz_tissue import visualize
with open(parametersFile) as jFile:
parameters = json.load(jFile)
visualize(parameters)
time.sleep(2)
\ No newline at end of file
basic_model_only_density/tissue.py
......@@ -5,6 +5,10 @@ import os
import h5py
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider
from tempfile import TemporaryDirectory
import vedo as vd
from matplotlib.animation import FuncAnimation
from mpl_toolkits.axes_grid1 import make_axes_locatable
df.set_log_level(df.LogLevel.ERROR)
df.parameters['form_compiler']['optimize'] = True
......@@ -33,6 +37,7 @@ class Tissue(object):
self.rho_boundary = 0.0
else:
self.rho_boundary = self.average_rho
self.bc = df.DirichletBC(self.function_space.sub(2),
df.Constant(self.rho_boundary),
'on_boundary')
......@@ -77,6 +82,7 @@ class Tissue(object):
base_rho = 0.0
else:
base_rho = self.average_rho
# rho0 = df.interpolate(df.Constant(self.average_rho),self.function_space.sub(2).collapse())
rho0 = df.interpolate(df.Expression(
'base_rho + cos(PI*x[0]/L)',
L=self.system_size,
......@@ -87,6 +93,7 @@ class Tissue(object):
# add noise
noise_rho = (self.noise_level
* (2 * np.random.random(rho0.vector().size()) - 1))
rho0.vector()[:] = rho0.vector()[:] + noise_rho
VFS = self.function_space.sub(1).collapse()
......@@ -158,88 +165,6 @@ class Tissue(object):
self.vFile.close()
self.rhoFile.close()
def viz_tissue(params, DIR=''):
savesteps = int(params['maxtime']/params['savetime'])
times = np.arange(savesteps+1) * params['savetime']
# Read mesh geometry from h5 file
var = 'density'
h5 = h5py.File(os.path.join(DIR,
'%s_%s.h5' % (params['timestamp'], var)), "r")
# should be in the loop if remeshing
topology = np.array(h5['%s/%s_0/mesh/topology'%(var,var)])
geometry = []
for i in range(len(times)):
geometry.append(np.array(h5['%s/%s_%d/mesh/geometry'%(var,var,i)]))
h5.close()
geometry = np.array(geometry)
geometry, zeros= np.dsplit(geometry, 2)
mesh = df.IntervalMesh(params['resolution'],
-params['system_size']/2,
params['system_size']/2)
# Read data
u = np.zeros((len(times), mesh.num_vertices(), 1))
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']))
# Reading data
print('Reading data...')
for steps in progressbar.progressbar(range(savesteps+1)):
mesh.coordinates()[:] = geometry[steps]
SFS = df.FunctionSpace(mesh, 'P', 1)
VFS = df.VectorFunctionSpace(mesh, 'P', 1)
ui, vi, rhoi = df.Function(VFS), df.Function(VFS), df.Function(SFS)
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(1, int(u_vec.shape[0])).T
v_vec = vi.compute_vertex_values(mesh)
v[steps] = v_vec.reshape(1, int(v_vec.shape[0])).T
rho[steps] = rhoi.compute_vertex_values(mesh)
uFile.close()
vFile.close()
rhoFile.close()
fig, axes = plt.subplots(2, 1, sharex=True, figsize=(8,8))
axes[-1].set_xlabel(r'$x$')
axes[0].set_ylabel(r'$\rho$')
axes[1].set_ylabel(r'$v$')
axes[0].set_xlim(np.min(geometry), np.max(geometry))
axes[0].set_ylim(np.min(rho), np.max(rho))
axes[1].set_ylim(np.min(v), np.max(v))
rhoplot, = axes[0].plot(geometry[0], rho[0], 'g-', ms=3)
velplot, = axes[1].plot(geometry[0], v[0], 'r-', ms=3)
def update(value):
ti = np.abs(times-value).argmin()
rhoplot.set_ydata(rho[ti])
rhoplot.set_xdata(geometry[ti])
velplot.set_ydata(v[ti])
velplot.set_xdata(geometry[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()
if __name__ == '__main__':
import json, datetime
......@@ -250,7 +175,7 @@ if __name__ == '__main__':
with open('parameters.json') as jsonFile:
params = json.load(jsonFile)
timestamp = '123456' #datetime.datetime.now().strftime("%d%m%y-%H%M%S")
timestamp = datetime.datetime.now().strftime("%d%m%y-%H%M%S")
params['timestamp'] = timestamp
tissue = Tissue(params)
......@@ -259,4 +184,3 @@ if __name__ == '__main__':
with open(params['timestamp'] + '_parameters.json', "w") as fp:
json.dump(params, fp, indent=4)
viz_tissue(params)
\ No newline at end of file
basic_model_only_density/viz_tissue.py 0 → 100644
import dolfin as df
import numpy as np
import vedo as vd
import os
import h5py
from matplotlib.widgets import Slider
import progressbar
import matplotlib.pyplot as plt
from tempfile import TemporaryDirectory
def visualize(params, DIR=''):
savesteps = int(params['maxtime']/params['savetime'])
times = np.arange(savesteps+1) * params['savetime']
# Read mesh geometry from h5 file
var = 'density'
h5 = h5py.File(os.path.join(DIR,
'%s_%s.h5' % (params['timestamp'], var)), "r")
# should be in the loop if remeshing
topology = np.array(h5['%s/%s_0/mesh/topology'%(var,var)])
geometry = []
for i in range(len(times)):
geometry.append(np.array(h5['%s/%s_%d/mesh/geometry'%(var,var,i)]))
h5.close()
geometry = np.array(geometry)
geometry, zeros= np.dsplit(geometry, 2)
mesh = df.IntervalMesh(params['resolution'],
- params['system_size']/2,
params['system_size']/2)
# Read data
u = np.zeros((len(times), mesh.num_vertices(), 1))
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']))
# Reading data
print('Reading data...')
for steps in progressbar.progressbar(range(savesteps+1)):
mesh.coordinates()[:] = geometry[steps]
SFS = df.FunctionSpace(mesh, 'P', 1)
VFS = df.VectorFunctionSpace(mesh, 'P', 1)
ui, vi, rhoi = df.Function(VFS), df.Function(VFS), df.Function(SFS)
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(1, int(u_vec.shape[0])).T
v_vec = vi.compute_vertex_values(mesh)
v[steps] = v_vec.reshape(1, int(v_vec.shape[0])).T
rho[steps] = rhoi.compute_vertex_values(mesh)
uFile.close()
vFile.close()
rhoFile.close()
# interactive plot
fig, axes = plt.subplots(2, 1, sharex=True, figsize=(8,8))
axes[-1].set_xlabel(r'$x$')
axes[0].set_ylabel(r'$\rho/\rho_0$')
axes[1].set_ylabel(r'$v/v_0$')
plt.rc('font', size=14)
axes[0].set_xlim(np.min(geometry), np.max(geometry))
axes[0].set_ylim(np.min(rho), np.max(rho))
axes[1].set_ylim(np.min(v), np.max(v))
rhoplot, = axes[0].plot(geometry[0], rho[0], 'g-', ms=3)
velplot, = axes[1].plot(geometry[0], v[0], 'r-', ms=3)
def update(value):
ti = np.abs(times-value).argmin()
rhoplot.set_ydata(rho[ti])
rhoplot.set_xdata(geometry[ti])
velplot.set_ydata(v[ti])
velplot.set_xdata(geometry[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()
# make movie
print('Saving movie-...')
FPS = 50
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))):
slider.set_val(times[tt])
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()
# L(t) vs t
length = geometry[:,-1,0]-geometry[:,0,0]
fig1, ax1 = plt.subplots(1, 1, figsize=(8, 8))
ax1.set_xlabel('$t$')
ax1.set_ylabel('$L(t)$')
ax1.set_xlim(np.min(times), np.max(times))
ax1.set_ylim(np.min(length), np.max(length)+1)
ax1.plot(times, length)
plt.savefig("%s.png" %params['timestamp'])
# plt.show()
np.save('%s_length.npy' %params['timestamp'], length)
if __name__ == '__main__':
import argparse, json
parser = argparse.ArgumentParser()
parser.add_argument('-j','--jsonfile', help='parameters file', required=True)
parser.add_argument('--onscreen', action=argparse.BooleanOptionalAction)
args = parser.parse_args()
with open(args.jsonfile) as jsonFile:
params = json.load(jsonFile)
visualize(params, DIR=os.path.dirname(args.jsonfile))
\ No newline at end of file
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