Skip to content

Jigyasa Watwani / growth-pattern-control

Go to a project
Commit df925e32 by Vijay Kumar Krishnamurthy
Options

First attempt at combined model with/without c

parent 14f0bc3d
Inline Side-by-side
Showing with 111 additions and 35 deletions
fixed_boundaries.py
...@@ -14,13 +14,20 @@ class FixedBoundaries(object): ...@@ -14,13 +14,20 @@ class FixedBoundaries(object):
# create mesh, mixed element and function space # create mesh, mixed element and function space
self.mesh = df.IntervalMesh(self.resolution, 0.0, 2.0*np.pi*self.system_size_by_2PI) self.mesh = df.IntervalMesh(self.resolution, 0.0, 2.0*np.pi*self.system_size_by_2PI)
scalar_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1) scalar_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1)
mixed_element = df.MixedElement([scalar_element, scalar_element, scalar_element]) # u, v, rho if self.morphogen:
# u, v, rho, c
mixed_element = df.MixedElement([scalar_element, scalar_element,
scalar_element, scalar_element])
else:
# u, v, rho
mixed_element = df.MixedElement([scalar_element, scalar_element, scalar_element])
# define function space with this mixed element # define function space with this mixed element
self.function_space = df.FunctionSpace(self.mesh, mixed_element) self.function_space = df.FunctionSpace(self.mesh, mixed_element)
# dirichlet boundaries for u, v # dirichlet boundaries for u, v
bc1 = df.DirichletBC(self.function_space.sub(0), df.Constant(0.0), 'on_boundary') # u,v at boundary = 0 # u,v at boundary = 0
bc1 = df.DirichletBC(self.function_space.sub(0), df.Constant(0.0), 'on_boundary')
bc2 = df.DirichletBC(self.function_space.sub(1), df.Constant(0.0), 'on_boundary') bc2 = df.DirichletBC(self.function_space.sub(1), df.Constant(0.0), 'on_boundary')
self.bc = ([bc1,bc2]) self.bc = ([bc1,bc2])
...@@ -32,40 +39,77 @@ class FixedBoundaries(object): ...@@ -32,40 +39,77 @@ class FixedBoundaries(object):
def advection(self, conc, vel, tconc): def advection(self, conc, vel, tconc):
return (df.inner((vel * conc).dx(0), tconc)) return (df.inner((vel * conc).dx(0), tconc))
def active_stress(self, rho): def active_stress(self, rho, c):
return self.lamda * rho/(rho + self.saturation_rho) return self.lamda * rho/(rho + self.saturation_rho) * c
def passive_stress(self, u, v): def passive_stress(self, u, v):
return (self.elasticity * u.dx(0) + self.viscosity * v.dx(0)) return (self.elasticity * u.dx(0) + self.viscosity * v.dx(0))
def stress(self, u, v, rho): def stress(self, u, v, rho, c):
return (self.passive_stress(u, v) + self.active_stress(rho)) return (self.passive_stress(u, v) + self.active_stress(rho, c))
def reaction_rho(self, rho, trho): def reaction_rho(self, rho, trho):
return self.turnover_rho * df.inner(rho - self.average_rho, trho) return self.turnover_rho * df.inner(rho - self.average_rho, trho)
def reaction_rho(self, c, tc):
return self.turnover_c * df.inner(c - self.average_c, tc)
def diffusion_reaction_rho(self, rho, trho): def diffusion_reaction_rho(self, rho, trho):
return (self.diffusion_rho * df.inner(rho.dx(0), trho.dx(0)) + self.reaction_rho(rho, trho)) return (self.diffusion_rho * df.inner(rho.dx(0), trho.dx(0)) + self.reaction_rho(rho, trho))
def setup_initial_conditions(self, u0, rho0): def diffusion_reaction_c(self, c, tc):
return (self.diffusion_c * df.inner(c.dx(0), tc.dx(0)) + self.reaction_rho(c, tc))
def setup_initial_conditions(self, u0=None, rho0=None, c0=None):
if u0==None and rho0==None and c0==None:
u0 = df.interpolate(df.Constant(0), self.function_space.sub(0).collapse())
rho0 = df.interpolate(df.Constant(self.average_rho), self.function_space.sub(2).collapse())
# add noise
noise_u = self.noise_level * (2*np.random.random(u0.vector().size())-1)
u0.vector()[:] = u0.vector()[:] + noise_u
noise_rho = self.noise_level * (2*np.random.random(rho0.vector().size())-1)
rho0.vector()[:] = rho0.vector()[:] + noise_rho
if self.morphogen:
c0 = df.interpolate(df.Constant(self.average_c), self.function_space.sub(3).collapse())
# add noise
noise_c = self.noise_level * (2*np.random.random(c0.vector().size())-1)
c0.vector()[:] = c0.vector()[:] + noise_c
else:
c0 = df.Constant(1.0)
else:
u0 = df.interpolate(u0, self.function_space.sub(0).collapse()) u0 = df.interpolate(u0, self.function_space.sub(0).collapse())
rho0 = df.interpolate(rho0, self.function_space.sub(2).collapse()) rho0 = df.interpolate(rho0, self.function_space.sub(2).collapse())
if self.morphogen:
c0 = df.interpolate(c0, self.function_space.sub(3).collapse())
else:
c0 = df.Constant(1.0)
VFS = self.function_space.sub(1).collapse() VFS = self.function_space.sub(1).collapse()
v = df.Function(VFS) v0 = df.Function(VFS)
tv = df.TestFunction(VFS) tv = df.TestFunction(VFS)
vform = (self.friction * df.inner(v, tv) vform = (self.friction * df.inner(v0, tv)
+ df.inner(self.stress(u0, v, rho0), tv.dx(0)) + df.inner(self.stress(u0, v0, rho0, c0), tv.dx(0))
) * df.dx ) * df.dx
bc = df.DirichletBC(VFS, df.Constant(0.0), 'on_boundary') bc = df.DirichletBC(VFS, df.Constant(0.0), 'on_boundary')
df.solve(vform == 0, v, bc) df.solve(vform == 0, v0, bc)
df.assign(self.f0, [u0, v, rho0]) if self.morphogen:
df.assign(self.f0, [u0, v0, rho0, c0])
else:
df.assign(self.f0, [u0, v0, rho0])
self.f1.assign(self.f0) self.f1.assign(self.f0)
def setup_weak_forms(self): def setup_weak_forms(self):
if self.morphogen:
u1, v1, rho1, c1 = df.split(self.f1)
u0, v0, rho0, c0 = df.split(self.f0)
u, v, rho, c = df.split(self.f)
tu, tv, trho, tc = df.TestFunctions(self.function_space)
else:
u1, v1, rho1 = df.split(self.f1) u1, v1, rho1 = df.split(self.f1)
u0, v0, rho0 = df.split(self.f0) u0, v0, rho0 = df.split(self.f0)
u, v, rho = df.split(self.f) u, v, rho = df.split(self.f)
tu, tv, trho = df.TestFunctions(self.function_space) tu, tv, trho = df.TestFunctions(self.function_space)
c = df.Constant(1.0)
uform = (df.inner((u - u0)/self.timestep, tu) uform = (df.inner((u - u0)/self.timestep, tu)
- 9/16 * df.inner(v, tu) - 9/16 * df.inner(v, tu)
...@@ -74,7 +118,7 @@ class FixedBoundaries(object): ...@@ -74,7 +118,7 @@ class FixedBoundaries(object):
) )
vform = (self.friction * df.inner(v, tv) vform = (self.friction * df.inner(v, tv)
+ df.inner(self.stress(u, v, rho), tv.dx(0)) + df.inner(self.stress(u, v, rho, c), tv.dx(0))
) )
rhoform = (df.inner((rho - rho0)/self.timestep, trho) rhoform = (df.inner((rho - rho0)/self.timestep, trho)
...@@ -84,33 +128,53 @@ class FixedBoundaries(object): ...@@ -84,33 +128,53 @@ class FixedBoundaries(object):
+ 3/8 * self.diffusion_reaction_rho(rho0, trho) + 3/8 * self.diffusion_reaction_rho(rho0, trho)
+ 1/16 * self.diffusion_reaction_rho(rho1, trho) + 1/16 * self.diffusion_reaction_rho(rho1, trho)
) )
self.form = (uform + vform + rhoform) * df.dx
def solve(self, u0, rho0): if self.morphogen:
times = np.arange(0.0, self.maxtime, self.timestep) cform = (df.inner((c - c0)/self.timestep, tc)
x = self.mesh.coordinates() + 3/2 * self.advection(c0, v, tc)
u_array = np.zeros((len(times), len(x))) - 1/2 * self.advection(c1, v, tc)
v_array = np.zeros_like(u_array) + 9/16 * self.diffusion_reaction_c(c, tc)
rho_array = np.zeros_like(u_array) + 3/8 * self.diffusion_reaction_c(c0, tc)
+ 1/16 * self.diffusion_reaction_c(c1, tc)
)
self.form = (uform + vform + rhoform + cform) * df.dx
else:
self.form = (uform + vform + rhoform) * df.dx
self.setup_initial_conditions(u0, rho0) 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 solve(self, extend=None, DIR=''):
if extend:
# read the last time point of earlier simulation
# u0, rho0, c0
pass
else:
u0, rho0, c0 = None, None, None
self.setup_initial_conditions(u0, rho0, c0)
self.setup_weak_forms() self.setup_weak_forms()
u, v, rho = self.f0.split(deepcopy=True) # store
u_array[0] = u.compute_vertex_values(self.mesh) data = []
v_array[0] = v.compute_vertex_values(self.mesh) data.append(self.get_data(self.f0))
rho_array[0] = rho.compute_vertex_values(self.mesh)
times = np.arange(0.0, self.maxtime, self.timestep)
for i in progressbar.progressbar(range(0, len(times))): for i in progressbar.progressbar(range(0, len(times))):
df.solve(self.form == 0, self.f, self.bc) df.solve(self.form == 0, self.f, self.bc)
u, v, rho = self.f.split(deepcopy=True) data.append(self.get_data(self.f))
u_array[i] = u.compute_vertex_values(self.mesh)
v_array[i] = v.compute_vertex_values(self.mesh)
rho_array[i] = rho.compute_vertex_values(self.mesh)
self.f0.assign(self.f) self.f0.assign(self.f)
self.f1.assign(self.f0) self.f1.assign(self.f0)
return (u_array, v_array, rho_array)
# parse data
data = np.array(data)
print(data.shape)
if self.morphogen:
# u, v, rho, c
return np.split(data, 4, axis=2)
else:
# u, v, rho
return np.split(data, 3, axis=2)
if __name__ == '__main__': if __name__ == '__main__':
...@@ -118,19 +182,27 @@ if __name__ == '__main__': ...@@ -118,19 +182,27 @@ if __name__ == '__main__':
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from matplotlib.widgets import Slider from matplotlib.widgets import Slider
with open('growth_parameters.json') as jsonFile: with open('fixed_boundaries.json') as jsonFile:
params = json.load(jsonFile) params = json.load(jsonFile)
fb = FixedBoundaries(params) fb = FixedBoundaries(params)
Z = fb.solve()
u0 = df.Expression('sin(x[0])', degree=1) if params['morphogen']:
rho0 = df.Expression('rho0 + 0.1 * (cos(x[0])+cos(x[0]/2))', rho0=fb.average_rho, degree=1) u, v, rho, c = Z
else:
u, v, rho = Z
(u, v, rho) = fb.solve(u0, rho0) print(u.shape, v.shape, rho.shape)
x = fb.mesh.coordinates() x = fb.mesh.coordinates()
times = np.arange(0, fb.maxtime, fb.timestep) times = np.arange(0, fb.maxtime, fb.timestep)
if params['morphogen']:
fig, (axu, axv, axrho, axc) = plt.subplots(4,1, sharex=True, figsize=(8,8))
axc.set_xlabel(r'$x$')
axc.set_ylim(np.min(c), np.max(c))
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)$')
...@@ -143,12 +215,16 @@ if __name__ == '__main__': ...@@ -143,12 +215,16 @@ if __name__ == '__main__':
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']:
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])
rhoplot.set_ydata(rho[ti]) rhoplot.set_ydata(rho[ti])
if params['morphogen']:
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])
......
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