Skip to content

Jigyasa Watwani / growth-pattern-control

Go to a project
Commit d7f785ec by Jigyasa Watwani
Options

working extend code

parent b4081442
Inline Side-by-side
Showing with 177 additions and 90 deletions
basic_model_only_density/parameters.json
......@@ -5,15 +5,15 @@
"viscosity": 1,
"zero_rho_boundary": true,
"friction": 1.0,
"lamda": -0.1,
"lamda": 1,
"diffusion_rho": 0.1,
"turnover_rho": 10.0,
"turnover_rho": 1,
"active_death": false,
"active_stress_setpoint": 1,
"active_stress_setpoint": 3,
"average_rho": 1.0,
"saturation_rho": 1.0,
"noise_level": 0.0,
"timestep": 0.001,
"timestep": 0.01,
"savetime": 0.1,
"maxtime": 70.0
"maxtime": 1.0
}
\ No newline at end of file
basic_model_only_density/run_tissue.py
......@@ -10,8 +10,7 @@ import numpy as np
parser = argparse.ArgumentParser()
parser.add_argument('-j','--jsonfile', help='data file',
default='parameters.json')
parser.add_argument('-t','--time', help='time to run', type=float,
default=100)
parser.add_argument('-t','--time', help='time to run', type=float)
args = parser.parse_args()
assert os.path.isfile(args.jsonfile), '%s file not found' % args.jsonfile
......@@ -19,68 +18,61 @@ with open(args.jsonfile) as jsonFile:
parameters = json.load(jsonFile)
if args.jsonfile=='parameters.json':
extend=False
extend = False
print('Fresh run...')
timestamp = datetime.datetime.now().strftime("%d%m%y-%H%M%S")
parameters['timestamp'] = timestamp
parametersFile = timestamp+'_parameters.json'
initu=None
initrho=None
oldMaxTime=0.0
parametersFile = timestamp + '_parameters.json'
initu = None
initrho = None
initTime = 0.0
mesh = None
else:
# extend=True
# print('Extending run %s...' % parameters['timestamp'])
# parametersFile = args.jsonfile
# # read geometry from h5 file
# var = 'density'
# h5 = h5py.File( '%s_%s.h5' % (parameters['timestamp'], var))
# # should be in the loop if remeshing
# savesteps = int(parameters['maxtime']/parameters['savetime'])
# geometry = np.array(h5['%s/%s_%d/mesh/geometry'%(var,var,savesteps)])
# # print(geometry[:,0])
# # geometry, zeros= np.dsplit(geometry, 2)
# mesh = df.IntervalMesh(parameters['resolution'],
# - parameters['system_size']/2,
# parameters['system_size']/2)
# print(mesh.coordinates())
# # Read data
# uFile = df.XDMFFile(
# '%s_displacement.xdmf' % parameters['timestamp'])
# rhoFile = df.XDMFFile(
# '%s_density.xdmf' % parameters['timestamp'])
# # Reading data
# print('Reading data...')
# mesh.coordinates()[:,0] = geometry[:,0]
# print(mesh.coordinates())
# SFS = df.FunctionSpace(mesh, 'P', 1)
# VFS = df.VectorFunctionSpace(mesh, 'P', 1)
# initu, initrho = df.Function(VFS), df.Function(SFS)
# uFile.read_checkpoint(initu, 'displacement', savesteps)
# rhoFile.read_checkpoint(initrho, 'density', savesteps)
# uFile.close()
# rhoFile.close()
# oldMaxTime = parameters['maxtime']
tissue = Tissue(parameters)
tissue.solve(initu, initrho, oldMaxTime)
extend = True
print('Extending run %s...' % parameters['timestamp'])
parametersFile = args.jsonfile
savesteps = int(parameters['maxtime']/parameters['savetime'])
mesh = df.IntervalMesh(parameters['resolution'],
- parameters['system_size']/2,
parameters['system_size']/2)
print('Before reading geometry the mesh:', mesh.coordinates()[:,0])
# read geometry from h5 file
var = 'density'
h5 = h5py.File( '%s_%s.h5' % (parameters['timestamp'], var,), 'a')
# print("Name of the file: ", h5.filename)
# print('Details of the h5 file:', [(i,j) for i,j in h5.attrs.items()])
print(f'var {var}, savesteps {savesteps}')
geometry = np.array(h5['%s/%s_%d/mesh/geometry'%(var,var, savesteps)])
mesh.coordinates()[:,0] = geometry[:,0]
print('After reading geometry the mesh:',mesh.coordinates()[:,0])
# Read data
SFS = df.FunctionSpace(mesh, 'P', 1)
VFS = df.VectorFunctionSpace(mesh, 'P', 1)
initu, initrho = df.Function(VFS), df.Function(SFS)
uFile = df.XDMFFile('%s_displacement.xdmf' % parameters['timestamp'])
rhoFile = df.XDMFFile('%s_density.xdmf' % parameters['timestamp'])
uFile.read_checkpoint(initu, 'displacement', savesteps)
rhoFile.read_checkpoint(initrho, 'density', savesteps)
uFile.close()
rhoFile.close()
initTime = parameters['maxtime']
print('old max time=', initTime)
parameters['maxtime'] = args.time
tissue = Tissue(parameters, mesh)
tissue.solve(initu, initrho, initTime, extend)
if extend:
parameters['maxtime'] = oldMaxTime + parameters['maxtime']
parameters['maxtime'] = initTime + parameters['maxtime']
with open(parametersFile, "w") as jsonFile:
json.dump(parameters, jsonFile, indent=4, sort_keys=True)
from viz_tissue import visualize
visualize(parameters, DIR='')
\ No newline at end of file
# from viz_tissue import visualize
# visualize(parameters, DIR='')
\ No newline at end of file
basic_model_only_density/tissue.py
import numpy as np
import dolfin as df
import progressbar
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
class Tissue(object):
def __init__(self, parameters):
def __init__(self, parameters, mesh=None):
# read in parameters
for key in parameters:
setattr(self, key, parameters[key])
self.mesh = df.Mesh()
# if no mesh is given
if mesh is None:
self.mesh = df.IntervalMesh(self.resolution,
-self.system_size/2,
self.system_size/2)
else:
self.mesh = mesh
print('Mesh passed as initial condition:', self.mesh.coordinates()[:,0])
scalar_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1)
vector_element = df.VectorElement('P', self.mesh.ufl_cell(), 1)
......@@ -71,10 +67,9 @@ class Tissue(object):
return (self.passive_stress(u, v) + self.active_stress(rho))
def diffusion_reaction_rho(self, rho, trho):
# NOTE: implement a step-function of density in the reaction
return (self.diffusion_rho * df.inner(df.nabla_grad(rho),
df.nabla_grad(trho))
- self.turnover_rho * df.inner(rho*(1 - rho/self.average_rho), # logistic growth, no need to implement step function
- self.turnover_rho * df.inner(rho*(1 - rho/self.average_rho),
trho)
)
......@@ -91,7 +86,6 @@ 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 + 0.1 * cos(PI*x[0]/L)*cos(PI*x[0]/L)',
L=self.system_size,
......@@ -138,7 +132,7 @@ class Tissue(object):
self.form = (uform + vform + rhoform) * df.dx
def solve(self, initu=None, initrho=None, inittime=None):
def solve(self, initu=None, initrho=None, initTime=0, extend=False):
self.uFile = df.XDMFFile(
'%s_displacement.xdmf' % self.timestamp)
self.vFile = df.XDMFFile(
......@@ -150,20 +144,26 @@ class Tissue(object):
self.setup_weak_forms()
# time-variables
self.time = inittime
self.time = initTime
savesteps = int(self.savetime/self.timestep)
maxsteps = int(self.maxtime/self.timestep)
# if this is a fresh run, write the initial condition
if not extend:
u, v, rho = self.function0.split(deepcopy=True)
self.uFile.write_checkpoint(u, 'displacement', self.time)
self.vFile.write_checkpoint(v, 'velocity', self.time)
self.rhoFile.write_checkpoint(rho, 'density', self.time)
for steps in progressbar.progressbar(range(1, maxsteps + 1)):
print(self.time)
df.solve(self.form == 0, self.function, self.bc)
self.function0.assign(self.function)
self.time += self.timestep
u, v, rho = self.function0.split(deepcopy=True)
if steps % savesteps == 0:
# if not extend:
self.uFile.write_checkpoint(u, 'displacement',
self.time, append=True)
self.vFile.write_checkpoint(v, 'velocity',
......@@ -174,7 +174,13 @@ class Tissue(object):
dr = df.project(v*self.timestep, self.function_space.sub(0).collapse())
df.ALE.move(self.mesh, dr)
self.uFile.close()
self.vFile.close()
self.rhoFile.close()
print('After solving till new maxtime, the mesh coordinates:',self.mesh.coordinates()[:,0])
self.uFile = None
self.vFile = None
self.rhoFile = None
# self.uFile.closed
# self.vFile.closed
# self.rhoFile.closed
basic_model_only_density/viz_tissue.py
import dolfin as df
import numpy as np
import vedo as vd
import os
import h5py
from matplotlib.widgets import Slider
......@@ -11,27 +10,33 @@ from tempfile import TemporaryDirectory
def visualize(params, DIR=''):
savesteps = int(params['maxtime']/params['savetime'])
times = np.arange(savesteps+1) * params['savetime']
print(times)
# 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)):
for i in progressbar.progressbar(range(len(times))):
print(times[i])
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)
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()))
elastic_stress = np.zeros((len(times), mesh.num_vertices()))
viscous_stress = np.zeros_like(elastic_stress)
active_stress = np.zeros_like(elastic_stress)
uFile = df.XDMFFile(os.path.join(DIR,
'%s_displacement.xdmf' % params['timestamp']))
......@@ -51,27 +56,47 @@ def visualize(params, DIR=''):
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)
e_stress = params['elasticity'] * ui.dx(0)
e_stress = df.project(e_stress, VFS)
elastic_stress[steps] = e_stress.compute_vertex_values(mesh)
v_stress = params['viscosity'] * vi.dx(0)
v_stress = df.project(v_stress, VFS)
viscous_stress[steps] = v_stress.compute_vertex_values(mesh)
a_stress = -params['lamda'] * rhoi/(rhoi + params['saturation_rho'])
a_stress = df.project(a_stress, SFS)
active_stress[steps] = a_stress.compute_vertex_values(mesh)
uFile.close()
vFile.close()
rhoFile.close()
# interactive plot
fig, axes = plt.subplots(2, 1, sharex=True, figsize=(8,8))
plt.rcParams.update({'font.size': 15})
fig, axes = plt.subplots(3, 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[2].set_ylabel(r'$u/u_0$')
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))
axes[0].set_ylim(np.min(rho)-0.1, np.max(rho)+0.1)
axes[1].set_ylim(np.min(v)-0.1, np.max(v)+0.1)
axes[2].set_ylim(np.min(u)-0.1, np.max(u)+0.1)
rhoplot, = axes[0].plot(geometry[0], rho[0], 'g-', ms=3)
velplot, = axes[1].plot(geometry[0], v[0], 'r-', ms=3)
uplot, = axes[2].plot(geometry[0], u[0], 'b-', ms=3)
def update(value):
ti = np.abs(times-value).argmin()
......@@ -79,6 +104,8 @@ def visualize(params, DIR=''):
rhoplot.set_xdata(geometry[ti])
velplot.set_ydata(v[ti])
velplot.set_xdata(geometry[ti])
uplot.set_ydata(u[ti])
uplot.set_xdata(geometry[ti])
plt.draw()
sax = plt.axes([0.1, 0.92, 0.7, 0.02])
......@@ -88,11 +115,63 @@ def visualize(params, DIR=''):
slider.drawon = False
slider.on_changed(update)
# plt.show()
print('Saving movie-...')
FPS = 50
movFile = os.path.join(DIR, '%s_fields.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()
# plotting the stresses
fig2, axes1 = plt.subplots(3, 1, sharex=True, figsize=(8,8))
axes1[-1].set_xlabel(r'$x$')
axes1[2].set_ylabel('Elastic stress')
axes1[1].set_ylabel('Viscous stress')
axes1[0].set_ylabel('Active stress')
axes1[0].set_xlim(np.min(geometry), np.max(geometry))
axes1[2].set_ylim(np.min(elastic_stress)-0.1, np.max(elastic_stress)+0.1)
axes1[1].set_ylim(np.min(viscous_stress)-0.1, np.max(viscous_stress)+0.1)
axes1[0].set_ylim(np.min(active_stress)-0.1, np.max(active_stress)+0.1)
e_stress_plot, = axes1[2].plot(geometry[0], elastic_stress[0], 'g-', ms=3)
v_stress_plot, = axes1[1].plot(geometry[0], viscous_stress[0], 'r-', ms=3)
a_stress_plot, = axes1[0].plot(geometry[0], -active_stress[0], 'b-', ms=3)
def update2(value):
ti = np.abs(times-value).argmin()
e_stress_plot.set_ydata(elastic_stress[ti])
e_stress_plot.set_xdata(geometry[ti])
v_stress_plot.set_ydata(viscous_stress[ti])
v_stress_plot.set_xdata(geometry[ti])
a_stress_plot.set_ydata(active_stress[ti])
a_stress_plot.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(update2)
# plt.show()
# make movie
print('Saving movie-...')
FPS = 50
movFile = os.path.join(DIR, '%s.mov' % params['timestamp'])
movFile = os.path.join(DIR, '%s_stresses.mov' % params['timestamp'])
fps = float(FPS)
command = "ffmpeg -y -r"
options = "-b:v 3600k -qscale:v 4 -vcodec mpeg4"
......@@ -101,7 +180,7 @@ def visualize(params, DIR=''):
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)
fig2.savefig(fr, facecolor=fig2.get_facecolor(), dpi=100)
os.system(command + " " + str(fps)
+ " -i " + tmp_dir.name + os.sep
+ "%03d.png " + options + " " + movFile)
......@@ -119,6 +198,15 @@ def visualize(params, DIR=''):
# plt.show()
np.save('%s_length.npy' %params['timestamp'], length)
fig2, ax2 = plt.subplots(1, 1, figsize=(8, 8))
ax2.set_xlabel('$t$')
ax2.set_ylabel('$L(t)$')
# ax2.set_xlim(np.min(times), np.max(times))
# ax2.set_ylim(np.min(length), np.max(length)+1)
ax2.loglog(times, length)
plt.savefig("%s_loglog.png" %params['timestamp'])
# plt.show()
if __name__ == '__main__':
import argparse, json
parser = argparse.ArgumentParser()
......
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