Skip to content

Jigyasa Watwani / growth-pattern-control

Go to a project
Commit d7dd851d by Jigyasa Watwani
Options

almost working extend code

parent 54022b0c
Inline Side-by-side
Showing with 98 additions and 62 deletions
basic_model_only_density/parameters.json
{ {
"resolution": 1000, "resolution": 3,
"system_size": 1.0, "system_size": 1.0,
"elasticity": 0.5, "elasticity": 1,
"viscosity": 1, "viscosity": 1,
"zero_rho_boundary": true, "zero_rho_boundary": true,
"friction": 1.0, "friction": 1.0,
"lamda": 1.5, "lamda": 1,
"diffusion_rho": 0.1, "diffusion_rho": 0.1,
"turnover_rho": 1, "turnover_rho": 1,
"active_death": false, "active_death": false,
...@@ -13,7 +13,8 @@ ...@@ -13,7 +13,8 @@
"average_rho": 1.0, "average_rho": 1.0,
"saturation_rho": 1.0, "saturation_rho": 1.0,
"noise_level": 0.0, "noise_level": 0.0,
"timestep": 0.01, "timestep": 0.1,
"savetime": 0.1, "savetime": 0.1,
"maxtime": 600.0 "maxtime": 0.5
} }
\ No newline at end of file
basic_model_only_density/run_tissue.py
...@@ -24,6 +24,7 @@ if args.jsonfile=='parameters.json': ...@@ -24,6 +24,7 @@ if args.jsonfile=='parameters.json':
parameters['timestamp'] = timestamp parameters['timestamp'] = timestamp
parametersFile = timestamp + '_parameters.json' parametersFile = timestamp + '_parameters.json'
initu = None initu = None
initv = None
initrho = None initrho = None
initTime = 0.0 initTime = 0.0
mesh = None mesh = None
...@@ -41,26 +42,28 @@ else: ...@@ -41,26 +42,28 @@ else:
# read geometry from h5 file # read geometry from h5 file
var = 'density' var = 'density'
h5 = h5py.File( '%s_%s.h5' % (parameters['timestamp'], var,), 'a') h5 = h5py.File( '%s_%s.h5' % (parameters['timestamp'], var,), 'r+')
geometry = np.array(h5['%s/%s_%d/mesh/geometry'%(var,var, savesteps)]) geometry = np.array(h5['%s/%s_%d/mesh/geometry'%(var,var, savesteps)])
mesh.coordinates()[:,0] = geometry[:,0] mesh.coordinates()[:,0] = geometry[:,0]
# Read data # Read data
SFS = df.FunctionSpace(mesh, 'P', 1) SFS = df.FunctionSpace(mesh, 'P', 1)
VFS = df.VectorFunctionSpace(mesh, 'P', 1) VFS = df.VectorFunctionSpace(mesh, 'P', 1)
initu, initrho = df.Function(VFS), df.Function(SFS) initu, initv, initrho = df.Function(VFS), df.Function(VFS), df.Function(SFS)
uFile = df.XDMFFile('%s_displacement.xdmf' % parameters['timestamp']) uFile = df.XDMFFile('%s_displacement.xdmf' % parameters['timestamp'])
vFile = df.XDMFFile('%s_velocity.xdmf' % parameters['timestamp'])
rhoFile = df.XDMFFile('%s_density.xdmf' % parameters['timestamp']) rhoFile = df.XDMFFile('%s_density.xdmf' % parameters['timestamp'])
uFile.read_checkpoint(initu, 'displacement', savesteps) uFile.read_checkpoint(initu, 'displacement', savesteps)
vFile.read_checkpoint(initv, 'velocity', savesteps)
rhoFile.read_checkpoint(initrho, 'density', savesteps) rhoFile.read_checkpoint(initrho, 'density', savesteps)
uFile.close() uFile.close()
vFile.close()
rhoFile.close() rhoFile.close()
initTime = parameters['maxtime'] initTime = parameters['maxtime']
parameters['maxtime'] = args.time parameters['maxtime'] = args.time
tissue = Tissue(parameters, mesh) tissue = Tissue(parameters, mesh)
tissue.solve(initu, initrho, initTime, extend) tissue.solve(initu, initv, initrho, initTime, extend)
if extend: if extend:
parameters['maxtime'] = initTime + parameters['maxtime'] parameters['maxtime'] = initTime + parameters['maxtime']
......
basic_model_only_density/tissue.py
...@@ -28,6 +28,7 @@ class Tissue(object): ...@@ -28,6 +28,7 @@ class Tissue(object):
# 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)
# boundary condition on rho
if self.zero_rho_boundary: if self.zero_rho_boundary:
self.rho_boundary = 0.0 self.rho_boundary = 0.0
else: else:
...@@ -37,6 +38,7 @@ class Tissue(object): ...@@ -37,6 +38,7 @@ class Tissue(object):
df.Constant(self.rho_boundary), df.Constant(self.rho_boundary),
'on_boundary') 'on_boundary')
# define functions on function space
self.function0 = df.Function(self.function_space) self.function0 = df.Function(self.function_space)
self.function = df.Function(self.function_space) self.function = df.Function(self.function_space)
...@@ -92,24 +94,19 @@ class Tissue(object): ...@@ -92,24 +94,19 @@ class Tissue(object):
degree=1, PI=np.pi), degree=1, PI=np.pi),
self.function_space.sub(2).collapse()) 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
else: else:
rho0 = df.interpolate(initrho, self.function_space.sub(2).collapse()) rho0 = df.interpolate(initrho, self.function_space.sub(2).collapse())
# velocity from u0 and rho0
VFS = self.function_space.sub(1).collapse() VFS = self.function_space.sub(1).collapse()
v0 = df.Function(VFS) v0 = df.Function(VFS)
tv = df.TestFunction(VFS) tv0 = df.TestFunction(VFS)
vform = (self.friction * df.inner(v0, tv) v0form = (self.friction * df.inner(v0, tv0)
+ df.inner(self.stress(u0, v0, rho0), + df.inner(self.stress(u0, v0, rho0),
self.epsilon(tv)) self.epsilon(tv0))
) * df.dx ) * df.dx
df.solve(vform == 0, v0) df.solve(v0form == 0, v0)
df.assign(self.function0, [u0, v0, rho0]) df.assign(self.function0, [u0, v0, rho0])
def setup_weak_forms(self): def setup_weak_forms(self):
...@@ -131,7 +128,8 @@ class Tissue(object): ...@@ -131,7 +128,8 @@ class Tissue(object):
self.form = (uform + vform + rhoform) * df.dx self.form = (uform + vform + rhoform) * df.dx
def solve(self, initu=None, initrho=None, initTime=0, extend=False): def solve(self, initu=None, initv=None, initrho=None, initTime=0, extend=False):
# create files if they don't exist, open them if they do exist
self.uFile = df.XDMFFile( self.uFile = df.XDMFFile(
'%s_displacement.xdmf' % self.timestamp) '%s_displacement.xdmf' % self.timestamp)
self.vFile = df.XDMFFile( self.vFile = df.XDMFFile(
...@@ -139,28 +137,54 @@ class Tissue(object): ...@@ -139,28 +137,54 @@ class Tissue(object):
self.rhoFile = df.XDMFFile( self.rhoFile = df.XDMFFile(
'%s_density.xdmf' % self.timestamp) '%s_density.xdmf' % self.timestamp)
self.setup_initial_conditions(initu, initrho) self.setup_initial_conditions(initu, initrho)
self.setup_weak_forms()
# time-variables # time-variables
self.time = initTime self.time = initTime
savesteps = int(self.savetime/self.timestep) savesteps = int(self.savetime/self.timestep)
maxsteps = int(self.maxtime/self.timestep) maxsteps = int(self.maxtime/self.timestep)
# if this is a fresh run, write the initial condition # if this is a fresh run, write the initial condition, then move the mesh using initial velocity
# otherwise just move the mesh using initial velocity
u0, v0, rho0 = self.function0.split(deepcopy=True)
print('time=', self.time)
print('u0[-1]=', u0.compute_vertex_values(self.mesh)[-1])
print('u0[0]=', u0.compute_vertex_values(self.mesh)[0])
print('u0[-2]=', u0.compute_vertex_values(self.mesh)[-2])
print('v0[-1]=', v0.compute_vertex_values(self.mesh)[-1])
print('rho0[-2]=', rho0.compute_vertex_values(self.mesh)[-2])
print('rho0[-1]=', rho0.compute_vertex_values(self.mesh)[-1])
print('rho0[0]=', rho0.compute_vertex_values(self.mesh)[0])
if not extend: if not extend:
u, v, rho = self.function0.split(deepcopy=True) self.uFile.write_checkpoint(u0, 'displacement', self.time)
self.uFile.write_checkpoint(u, 'displacement', self.time) self.vFile.write_checkpoint(v0, 'velocity', self.time)
self.vFile.write_checkpoint(v, 'velocity', self.time) self.rhoFile.write_checkpoint(rho0, 'density', self.time)
self.rhoFile.write_checkpoint(rho, 'density', self.time) print('time=', self.time)
print('u0[-1]=', u0.compute_vertex_values(self.mesh)[-1])
print('u0[0]=', u0.compute_vertex_values(self.mesh)[0])
print('u0[-2]=', u0.compute_vertex_values(self.mesh)[-2])
print('v0[-1]=', v0.compute_vertex_values(self.mesh)[-1])
print('rho0[-2]=', rho0.compute_vertex_values(self.mesh)[-1])
print('rho0[-1]=', rho0.compute_vertex_values(self.mesh)[-1])
print('rho0[0]=', rho0.compute_vertex_values(self.mesh)[0])
# move mesh at the first time-point
dr0 = df.project(v0*self.timestep, self.function_space.sub(0).collapse())
df.ALE.move(self.mesh, dr0)
self.setup_weak_forms()
for steps in progressbar.progressbar(range(1, maxsteps + 1)): for steps in progressbar.progressbar(range(1, maxsteps + 1)):
self.time += self.timestep
df.solve(self.form == 0, self.function, self.bc) df.solve(self.form == 0, self.function, self.bc)
self.function0.assign(self.function) self.function0.assign(self.function)
self.time += self.timestep
u, v, rho = self.function0.split(deepcopy=True) u, v, rho = self.function0.split(deepcopy=True)
if steps % savesteps == 0: if steps % savesteps == 0:
print('time=',self.time)
# if not extend: # if not extend:
self.uFile.write_checkpoint(u, 'displacement', self.uFile.write_checkpoint(u, 'displacement',
self.time, append=True) self.time, append=True)
...@@ -168,12 +192,20 @@ class Tissue(object): ...@@ -168,12 +192,20 @@ class Tissue(object):
self.time, append=True) self.time, append=True)
self.rhoFile.write_checkpoint(rho, 'density', self.rhoFile.write_checkpoint(rho, 'density',
self.time, append=True) self.time, append=True)
print('u[-1]=', u.compute_vertex_values(self.mesh)[-1])
print('u[0]=', u.compute_vertex_values(self.mesh)[0])
print('u[-2]=', u.compute_vertex_values(self.mesh)[-2])
print('v[-1]=', v.compute_vertex_values(self.mesh)[-1])
print('rho[-2]=', rho.compute_vertex_values(self.mesh)[-2])
print('rho[-1]=', rho.compute_vertex_values(self.mesh)[-1])
print('rho[0]=', rho.compute_vertex_values(self.mesh)[0])
# move mesh # move mesh
dr = df.project(v*self.timestep, self.function_space.sub(0).collapse()) dr = df.project(v*self.timestep, self.function_space.sub(0).collapse())
df.ALE.move(self.mesh, dr) df.ALE.move(self.mesh, dr)
self.uFile.close() self.uFile.close()
self.vFile.close() self.vFile.close()
self.rhoFile.close() self.rhoFile.close()
basic_model_only_density/viz_tissue.py
...@@ -113,22 +113,22 @@ def visualize(params, DIR=''): ...@@ -113,22 +113,22 @@ def visualize(params, DIR=''):
slider.drawon = False slider.drawon = False
slider.on_changed(update) slider.on_changed(update)
print('Saving movie-...') # print('Saving movie-...')
FPS = 50 # FPS = 50
movFile = os.path.join(DIR, '%s_fields.mov' % params['timestamp']) # movFile = os.path.join(DIR, '%s_fields.mov' % params['timestamp'])
fps = float(FPS) # fps = float(FPS)
command = "ffmpeg -y -r" # command = "ffmpeg -y -r"
options = "-b:v 3600k -qscale:v 4 -vcodec mpeg4" # options = "-b:v 3600k -qscale:v 4 -vcodec mpeg4"
tmp_dir = TemporaryDirectory() # tmp_dir = TemporaryDirectory()
get_filename = lambda x: os.path.join(tmp_dir.name, x) # get_filename = lambda x: os.path.join(tmp_dir.name, x)
for tt in progressbar.progressbar(range(len(times))): # for tt in progressbar.progressbar(range(len(times))):
slider.set_val(times[tt]) # slider.set_val(times[tt])
fr = get_filename("%03d.png" % tt) # fr = get_filename("%03d.png" % tt)
fig.savefig(fr, facecolor=fig.get_facecolor(), dpi=100) # fig.savefig(fr, facecolor=fig.get_facecolor(), dpi=100)
os.system(command + " " + str(fps) # os.system(command + " " + str(fps)
+ " -i " + tmp_dir.name + os.sep # + " -i " + tmp_dir.name + os.sep
+ "%03d.png " + options + " " + movFile) # + "%03d.png " + options + " " + movFile)
tmp_dir.cleanup() # tmp_dir.cleanup()
# plotting the stresses # plotting the stresses
figure, axes1 = plt.subplots(3, 1, sharex=True, figsize=(8,8)) figure, axes1 = plt.subplots(3, 1, sharex=True, figsize=(8,8))
...@@ -167,26 +167,26 @@ def visualize(params, DIR=''): ...@@ -167,26 +167,26 @@ def visualize(params, DIR=''):
# plt.show() # plt.show()
# make movie # make movie
print('Saving movie-...') # print('Saving movie-...')
FPS = 50 # FPS = 50
movFile = os.path.join(DIR, '%s_stresses.mov' % params['timestamp']) # movFile = os.path.join(DIR, '%s_stresses.mov' % params['timestamp'])
fps = float(FPS) # fps = float(FPS)
command = "ffmpeg -y -r" # command = "ffmpeg -y -r"
options = "-b:v 3600k -qscale:v 4 -vcodec mpeg4" # options = "-b:v 3600k -qscale:v 4 -vcodec mpeg4"
tmp_dir = TemporaryDirectory() # tmp_dir = TemporaryDirectory()
get_filename = lambda x: os.path.join(tmp_dir.name, x) # get_filename = lambda x: os.path.join(tmp_dir.name, x)
for tt in progressbar.progressbar(range(len(times))): # for tt in progressbar.progressbar(range(len(times))):
slider.set_val(times[tt]) # slider.set_val(times[tt])
fr = get_filename("%03d.png" % tt) # fr = get_filename("%03d.png" % tt)
figure.savefig(fr, facecolor=figure.get_facecolor(), dpi=100) # figure.savefig(fr, facecolor=figure.get_facecolor(), dpi=100)
os.system(command + " " + str(fps) # os.system(command + " " + str(fps)
+ " -i " + tmp_dir.name + os.sep # + " -i " + tmp_dir.name + os.sep
+ "%03d.png " + options + " " + movFile) # + "%03d.png " + options + " " + movFile)
tmp_dir.cleanup() # tmp_dir.cleanup()
# L(t) vs t # L(t) vs t
length = geometry[:,-1,0]-geometry[:,0,0] length = geometry[:,-1,0]-geometry[:,0,0]
print(length) # print(length[-1])
fig1, ax1 = plt.subplots(1, 1, figsize=(8, 8)) fig1, ax1 = plt.subplots(1, 1, figsize=(8, 8))
ax1.set_xlabel('$t$') ax1.set_xlabel('$t$')
ax1.set_ylabel('$L(t)$') ax1.set_ylabel('$L(t)$')
...@@ -210,7 +210,6 @@ def visualize(params, DIR=''): ...@@ -210,7 +210,6 @@ def visualize(params, DIR=''):
# dldt vs t # dldt vs t
dldt = np.gradient(length, times) dldt = np.gradient(length, times)
print(dldt)
fig11, ax11 = plt.subplots(1, 1, figsize=(8, 8)) fig11, ax11 = plt.subplots(1, 1, figsize=(8, 8))
ax11.set_xlabel('$t$') ax11.set_xlabel('$t$')
ax11.set_ylabel('$dL/dt$') ax11.set_ylabel('$dL/dt$')
......
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