Skip to content

Jigyasa Watwani / growth-pattern-control

Go to a project
Commit 4a038a83 by Jigyasa Watwani
Options

passive stress field instead of total stress field, and d-dim derivatives

parent 920268c0
Inline Side-by-side
Showing with 180 additions and 155 deletions
basic_model_only_density/tissue_maxwell.py
......@@ -24,106 +24,141 @@ class Tissue(object):
self.system_size/2)
scalar_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1)
vector_element = df.VectorElement('P', self.mesh.ufl_cell(), 1)
tensor_element = df.TensorElement('P', self.mesh.ufl_cell(), 1)
# stress, v, rho
mixed_element = df.MixedElement([scalar_element, scalar_element, scalar_element])
# v, p_stress, rho
mixed_element = df.MixedElement([vector_element,
tensor_element,
scalar_element])
# define function space with this mixed element
self.function_space = df.FunctionSpace(self.mesh, mixed_element)
self.bc = df.DirichletBC(self.function_space.sub(2), df.Constant(0.0), 'on_boundary')
self.bc = df.DirichletBC(self.function_space.sub(2),
df.Constant(0.0),
'on_boundary'),
self.function0 = df.Function(self.function_space)
self.function = df.Function(self.function_space)
self.relaxation_time = self.viscosity/self.elasticity
self.tau = self.viscosity/self.elasticity
def advection(self, conc, vel, tconc):
return df.inner(df.div(vel * conc), tconc)
def active_stress(self, rho):
if self.active_death:
return (-self.lamda * rho * (rho - self.active_stress_setpoint)/(rho**2 + self.saturation_rho**2)
* df.Identity(1))
else:
return (-self.lamda * rho /(rho + self.saturation_rho)
* df.Identity(1))
def epsilon(self, v):
return df.sym(df.nabla_grad(v))
def stress(self, pstress, rho):
return (pstress + self.active_stress(rho))
def diffusion_reaction_rho(self, rho, trho):
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), trho)
)
def setup_initial_conditions(self):
# ic for stress
rho0 = df.interpolate(df.Expression('cos(pi * x[0]/L) * cos(pi * x[0]/L)', pi = np.pi, L = self.system_size, degree=1), self.function_space.sub(2).collapse())
stress0 = df.interpolate(df.Constant(0.0), self.function_space.sub(0).collapse())
SFS = self.function_space.sub(1).collapse()
v0 = df.Function(SFS)
tv0 = df.TestFunction(SFS)
zero_tensor = df.Constant(((0.0,),))
pstress0 = df.interpolate(zero_tensor,
self.function_space.sub(1).collapse())
rho0 = df.interpolate(df.Expression(
'0.1 * cos(PI*x[0]/L)*cos(PI*x[0]/L)',
L=self.system_size,
degree=1, PI=np.pi),
self.function_space.sub(2).collapse())
# 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(0).collapse()
v0 = df.Function(VFS)
tv0 = df.TestFunction(VFS)
v0form = (self.friction * df.inner(v0, tv0)
+ df.inner(stress0, tv0.dx(0))
) * df.dx
+ df.inner(self.stress(pstress0, rho0),
self.epsilon(tv0))
) * df.dx
df.solve(v0form == 0, v0)
df.assign(self.function0, [stress0, v0, rho0])
def active_stress(self, rho):
return -self.lamda * rho/(rho + self.saturation_rho)
df.assign(self.function0, [v0, pstress0, rho0])
def setup_weak_forms(self):
stress0, v0, rho0 = df.split(self.function0)
stress, v, rho = df.split(self.function)
tstress, tv, trho = df.TestFunctions(self.function_space)
vform = self.friction * df.inner(v, tv) + df.inner(stress, tv.dx(0))
stressform = ( df.inner(stress - self.active_stress(rho), tstress)
+ self.relaxation_time *
df.inner((1/self.timestep) * (stress - stress0 - self.active_stress(rho) + self.active_stress(rho0)), tstress)
+ self.relaxation_time *
df.inner(v0 * (stress0 - self.active_stress(rho0)).dx(0), tstress)
- (self.viscosity) *
df.inner(v.dx(0), tstress)
)
v0, pstress0, rho0 = df.split(self.function0)
v, pstress, rho = df.split(self.function)
tv, tpstress, trho = df.TestFunctions(self.function_space)
vform = (self.friction * df.inner(v, tv)
+ df.inner(self.stress(pstress, rho),
self.epsilon(tv)))
rhoform = (df.inner((rho - rho0)/self.timestep, trho)
+ df.inner((rho0 * v0).dx(0), trho)
+ self.diffusion_rho * df.inner(rho.dx(0), trho.dx(0))
- self.turnover_rho * df.inner(rho * (1 - rho/self.average_rho), trho)
)
+ self.advection(rho0, v0, trho)
+ self.diffusion_reaction_rho(rho, trho))
pstressform = ( self.tau * df.inner((pstress - pstress0)/self.timestep, tpstress)
+ self.tau * df.inner(df.dot(v, df.nabla_grad(pstress)), tpstress)
+ df.inner(pstress, tpstress)
- self.viscosity * df.inner(self.epsilon(v), tpstress)
)
self.form = (vform + stressform + rhoform)*df.dx
self.form = (pstressform + vform + rhoform) * df.dx
def solve(self, DIR=''):
self.stressFile = df.XDMFFile(os.path.join(DIR,
'%s_stress.xdmf' % self.timestamp))
self.vFile = df.XDMFFile(os.path.join(DIR,
self.pstressFile = df.XDMFFile(os.path.join(DIR,
'%s_pstress.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))
self.setup_initial_conditions()
self.setup_weak_forms()
# time-variables
self.time = 0.0
savesteps = int(self.savetime/self.timestep)
maxsteps = int(self.maxtime/self.timestep)
stress, v, rho = self.function0.split(deepcopy=True)
self.stressFile.write_checkpoint(stress, 'stress', self.time)
self.rhoFile.write_checkpoint(rho, 'density', self.time)
v, pstress, rho = self.function0.split(deepcopy=True)
self.pstressFile.write_checkpoint(pstress, 'pstress', 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)):
df.solve(self.form == 0, self.function, self.bc)
self.function0.assign(self.function)
self.time += self.timestep
stress, v, rho = self.function0.split(deepcopy=True)
v, pstress, rho = self.function0.split(deepcopy=True)
if steps % savesteps == 0:
self.stressFile.write_checkpoint(stress, 'stress',
self.pstressFile.write_checkpoint(pstress, 'pstress',
self.time, append=True)
self.rhoFile.write_checkpoint(rho, 'density',
self.time, append=True)
self.vFile.write_checkpoint(v, 'velocity',
self.time, append=True)
self.rhoFile.write_checkpoint(rho, 'density',
self.time, append=True)
# move mesh
dr = df.project(v * self.timestep, self.function_space.sub(1).collapse())
dr = df.project(v*self.timestep, self.function_space.sub(0).collapse())
df.ALE.move(self.mesh, dr)
self.stressFile.close()
self.pstressFile.close()
self.vFile.close()
self.rhoFile.close()
if __name__ == '__main__':
import json, datetime
......@@ -141,8 +176,7 @@ if __name__ == '__main__':
tissue.solve()
with open(params['timestamp'] + '_parameters.json', "w") as fp:
json.dump(params, fp, indent=4)
json.dump(params, fp, indent=4)
from viz_tissue_maxwell import visualize
visualize(params, DIR='')
basic_model_only_density/viz_tissue_maxwell.py
......@@ -4,10 +4,9 @@ import vedo as vd
import os
import h5py
from matplotlib.widgets import Slider
import progressbar
# import progressbar
import matplotlib.pyplot as plt
from tempfile import TemporaryDirectory
import numdifftools as nd
def visualize(params, DIR=''):
savesteps = int(params['maxtime']/params['savetime'])
......@@ -24,57 +23,65 @@ def visualize(params, DIR=''):
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)
- params['system_size']/2,
params['system_size']/2)
# Read data
pstress = np.zeros((len(times), mesh.num_vertices(), 2))
v = np.zeros((len(times), mesh.num_vertices(), 1))
rho = np.zeros((len(times), mesh.num_vertices()))
stress = np.zeros_like(rho)
pstressFile = df.XDMFFile(os.path.join(DIR,
'%s_pstress.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']))
stressFile = df.XDMFFile(os.path.join(DIR,
'%s_stress.xdmf' % params['timestamp']))
# Reading data
print('Reading data...')
for steps in progressbar.progressbar(range(savesteps+1)):
for steps in range(savesteps+1):
mesh.coordinates()[:] = geometry[steps]
SFS = df.FunctionSpace(mesh, 'P', 1)
rhoi, stressi = df.Function(SFS), df.Function(SFS)
VFS = df.VectorFunctionSpace(mesh, 'P', 1)
TFS = df.TensorFunctionSpace(mesh, 'P', 1)
pstressi, vi, rhoi = df.Function(TFS), df.Function(VFS), df.Function(SFS)
pstressFile.read_checkpoint(pstressi, 'pstress', steps)
vFile.read_checkpoint(vi, 'velocity', steps)
rhoFile.read_checkpoint(rhoi, 'density', steps)
stressFile.read_checkpoint(stressi, 'stress', steps)
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)
stress[steps] = stressi.compute_vertex_values(mesh)
pstressFile.close()
vFile.close()
rhoFile.close()
stressFile.close()
print(geometry[-1])
# interactive plot
plt.rcParams.update({'font.size': 18})
plt.yticks(fontsize=18)
plt.xticks(fontsize=18)
plt.rcParams.update({'font.size': 15})
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('stress')
axes[1].set_ylabel(r'$v/v_0$')
axes[0].set_xlim(np.min(geometry), np.max(geometry))
axes[0].set_ylim(0, np.max(rho)+0.05)
axes[1].set_ylim(np.min(stress)-0.05, np.max(stress)+0.05)
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)
rhoplot, = axes[0].plot(geometry[0], rho[0], 'g-', ms=3)
velplot, = axes[1].plot(geometry[0], stress[0], 'r-', ms=3)
velplot, = axes[1].plot(geometry[0], v[0], 'r-', ms=3)
def update(value):
ti = np.abs(times - value).argmin()
ti = np.abs(times-value).argmin()
rhoplot.set_ydata(rho[ti])
rhoplot.set_xdata(geometry[ti])
velplot.set_ydata(stress[ti])
velplot.set_ydata(v[ti])
velplot.set_xdata(geometry[ti])
plt.draw()
......@@ -84,58 +91,61 @@ def visualize(params, DIR=''):
fc='#999999')
slider.drawon = False
slider.on_changed(update)
plt.show()
# print('Saving movie-...')
# FPS = 100
# 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()
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 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
active_stress = -params['lamda'] * rho/(rho + params['saturation_rho'])
passive_stress = stress - active_stress
fig1, axes1 = plt.subplots(2, 1, sharex=True, figsize=(8,8))
axes1[-1].set_xlabel(r'$x$')
axes1[1].set_ylabel('Passive stress')
axes1[0].set_ylabel('Active stress')
axes1[0].set_xlim(np.min(geometry), np.max(geometry))
axes1[1].set_ylim(np.min(passive_stress)-0.05, np.max(passive_stress)+0.05)
axes1[0].set_ylim(np.min(active_stress)-0.05, np.max(active_stress)+0.05)
# 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)
p_stress_plot, = axes1[1].plot(geometry[0], passive_stress[0], 'b-', ms=3)
a_stress_plot, = axes1[0].plot(geometry[0], active_stress[0], 'r-', ms=3)
# 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 update_2(value):
ti = np.abs(times-value).argmin()
# def update2(value):
# ti = np.abs(times-value).argmin()
p_stress_plot.set_ydata(passive_stress[ti])
p_stress_plot.set_xdata(geometry[ti])
# e_stress_plot.set_ydata(elastic_stress[ti])
# e_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(update_2)
plt.show()
# v_stress_plot.set_ydata(viscous_stress[ti])
# v_stress_plot.set_xdata(geometry[ti])
# make movie
# 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_stresses.mov' % params['timestamp'])
......@@ -144,46 +154,27 @@ def visualize(params, DIR=''):
# 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))):
# for tt in range(len(times)):
# slider.set_val(times[tt])
# fr = get_filename("%03d.png" % tt)
# fig1.savefig(fr, facecolor=fig1.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)
# tmp_dir.cleanup()
# L(t) vs t
# length = geometry[:,-1,0]-geometry[:,0,0]
# 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)-0.05, np.max(length)+0.05)
# ax2.plot(times, length)
# plt.savefig("%s.png" %params['timestamp'])
# # plt.show()
# np.save('%s_length.npy' %params['timestamp'], length)
# print(length[-1])
# print(length[-1]-length[-3])
# fig3, ax3 = plt.subplots(1, 1, figsize=(8, 8))
# ax3.set_xlabel('$\log t$')
# ax3.set_ylabel('$\log L(t)$')
# ax3.loglog(times, length,'o')
# plt.savefig("%s_logLlogt.png" %params['timestamp'])
# fig4, ax4 = plt.subplots(1, 1, figsize=(8, 8))
# ax4.set_xlabel('$t$')
# ax4.set_ylabel('$\log L(t)$')
# ax4.semilogy(times, length,'o')
# plt.savefig("%s_logLt.png" %params['timestamp'])
# fig5, ax5 = plt.subplots(1, 1, figsize=(8, 8))
# ax5.set_xlabel('$\log t$')
# ax5.set_ylabel('$L(t)$')
# ax5.semilogx(times, length,'o')
# plt.savefig("%s_logtL.png" %params['timestamp'])
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)
print(length[-1])
if __name__ == '__main__':
import argparse, json
......
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