Skip to content

Jigyasa Watwani / growth-pattern-control

Go to a project
Commit 2d1f87c1 by Jigyasa Watwani
Options

fixed small bugs

parent 43c1d1e3
Inline Side-by-side
Showing with 55 additions and 64 deletions
basic_model_only_density/parameters.json
{
"resolution": 3,
"resolution": 128,
"system_size": 1.0,
"elasticity": 1,
"viscosity": 1,
......@@ -9,12 +9,12 @@
"diffusion_rho": 0.1,
"turnover_rho": 1,
"active_death": false,
"active_stress_setpoint": 3,
"active_stress_setpoint": 1,
"average_rho": 1.0,
"saturation_rho": 1.0,
"noise_level": 0.0,
"timestep": 0.1,
"timestep": 0.01,
"savetime": 0.1,
"maxtime": 0.5
"maxtime": 100
}
\ No newline at end of file
basic_model_only_density/run_tissue.py
......@@ -34,7 +34,7 @@ else:
print('Extending run %s...' % parameters['timestamp'])
parametersFile = args.jsonfile
savesteps = int(parameters['maxtime']/parameters['savetime'])
savesteps = round(parameters['maxtime']/parameters['savetime'])
mesh = df.IntervalMesh(parameters['resolution'],
- parameters['system_size']/2,
......
basic_model_only_density/tissue.py
......@@ -10,7 +10,8 @@ class Tissue(object):
# read in parameters
for key in parameters:
setattr(self, key, parameters[key])
# if no mesh is given
# if no mesh is given as initial condition, create one
if mesh is None:
self.mesh = df.IntervalMesh(self.resolution,
-self.system_size/2,
......@@ -20,6 +21,7 @@ class Tissue(object):
scalar_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1)
vector_element = df.VectorElement('P', self.mesh.ufl_cell(), 1)
# u, v, rho
mixed_element = df.MixedElement([vector_element,
vector_element,
......@@ -74,15 +76,11 @@ class Tissue(object):
trho)
)
def setup_initial_conditions(self, initu=None, initrho=None):
if initu == None:
def setup_initial_conditions(self):
zero_vector = df.Constant((0.0,))
u0 = df.interpolate(zero_vector,
self.function_space.sub(0).collapse())
else:
u0 = df.interpolate(initu, self.function_space.sub(0).collapse())
if initrho == None:
if self.zero_rho_boundary:
base_rho = 0.0
else:
......@@ -94,8 +92,6 @@ class Tissue(object):
degree=1, PI=np.pi),
self.function_space.sub(2).collapse())
else:
rho0 = df.interpolate(initrho, self.function_space.sub(2).collapse())
# velocity from u0 and rho0
VFS = self.function_space.sub(1).collapse()
......@@ -129,6 +125,7 @@ class Tissue(object):
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(
'%s_displacement.xdmf' % self.timestamp)
......@@ -137,72 +134,48 @@ class Tissue(object):
self.rhoFile = df.XDMFFile(
'%s_density.xdmf' % self.timestamp)
self.setup_initial_conditions(initu, initrho)
# time-variables
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, 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])
savesteps = round(self.savetime/self.timestep)
maxsteps = round(self.maxtime/self.timestep)
# if not extend, find the initial velocity, save the ic, move the mesh with this initial velocity
if not extend:
self.setup_initial_conditions()
u0, v0, rho0 = self.function0.split(deepcopy=True)
self.uFile.write_checkpoint(u0, 'displacement', self.time)
self.vFile.write_checkpoint(v0, 'velocity', self.time)
self.rhoFile.write_checkpoint(rho0, '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())
dr0 = df.project(v0 * self.timestep, self.function_space.sub(0).collapse())
df.ALE.move(self.mesh, dr0)
# if extend, move the mesh with this initial velocity
else:
u0 = df.project(initu, self.function_space.sub(0).collapse())
v0 = df.project(initv, self.function_space.sub(1).collapse())
rho0 = df.project(initrho, self.function_space.sub(2).collapse())
dr0 = df.project(v0 * self.timestep, self.function_space.sub(0).collapse())
df.ALE.move(self.mesh, dr0)
df.assign(self.function0, [u0, v0, rho0])
self.setup_weak_forms()
for steps in progressbar.progressbar(range(1, maxsteps + 1)):
self.time += self.timestep
df.solve(self.form == 0, self.function, self.bc)
self.function0.assign(self.function)
u, v, rho = self.function0.split(deepcopy=True)
if steps % savesteps == 0:
print('time=',self.time)
# if not extend:
self.uFile.write_checkpoint(u, 'displacement',
self.time, append=True)
self.vFile.write_checkpoint(v, 'velocity',
self.time, append=True)
self.rhoFile.write_checkpoint(rho, 'density',
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
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)
......
basic_model_only_density/viz_tissue.py
......@@ -8,13 +8,12 @@ import matplotlib.pyplot as plt
from tempfile import TemporaryDirectory
def visualize(params, DIR=''):
savesteps = int(params['maxtime']/params['savetime'])
savesteps = round(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")
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)])
......@@ -112,6 +111,7 @@ def visualize(params, DIR=''):
fc='#999999')
slider.drawon = False
slider.on_changed(update)
plt.show()
# print('Saving movie-...')
# FPS = 50
......@@ -164,7 +164,7 @@ def visualize(params, DIR=''):
fc='#999999')
slider.drawon = False
slider.on_changed(update2)
# plt.show()
plt.show()
# make movie
# print('Saving movie-...')
......@@ -192,20 +192,28 @@ def visualize(params, DIR=''):
ax1.set_ylabel('$L(t)$')
ax1.set_xlim(np.min(times), np.max(times))
ax1.set_ylim(np.min(length), np.max(length)+0.01)
ax1.plot(times, length)
ax1.tick_params(axis='both', which='both', labelsize=10)
ax1.grid(True)
ax1.plot(times, length, marker='o', ms=3, linestyle='None')
plt.savefig("%s.png" %params['timestamp'])
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.loglog(times, length)
ax2.tick_params(axis='both', which='both', labelsize=10)
ax2.grid(True)
ax2.loglog(times, length, marker='o', ms=3, linestyle='None')
plt.savefig("%s_loglog.png" %params['timestamp'])
fig3, ax3 = plt.subplots(1, 1, figsize=(8, 8))
ax3.set_xlabel('$t$')
ax3.set_ylabel('$L(t)$')
ax3.semilogy(times, length)
ax3.tick_params(axis='both', which='both', labelsize=10)
ax3.grid(True)
ax3.semilogy(times, length, marker='o', ms=3, linestyle='None')
plt.savefig("%s_semilog.png" %params['timestamp'])
# dldt vs t
......@@ -215,19 +223,28 @@ def visualize(params, DIR=''):
ax11.set_ylabel('$dL/dt$')
ax11.set_xlim(np.min(times), np.max(times))
ax11.set_ylim(np.min(dldt), np.max(dldt)+0.01)
ax11.plot(times, dldt)
ax11.tick_params(axis='both', which='both', labelsize=10)
ax11.grid(True)
ax11.plot(times, dldt, marker='o', ms=3, linestyle='None')
plt.savefig("%s_dldt.png" %params['timestamp'])
fig21, ax21 = plt.subplots(1, 1, figsize=(8, 8))
ax21.set_xlabel('$t$')
ax21.set_ylabel('$dL/dt$')
ax21.loglog(times, dldt)
ax21.tick_params(axis='both', which='both', labelsize=10)
ax21.grid(True)
ax21.loglog(times, dldt, marker='o', ms=3, linestyle='None')
plt.savefig("%s_loglog_dldt.png" %params['timestamp'])
fig31, ax31 = plt.subplots(1, 1, figsize=(8, 8))
ax31.set_xlabel('$t$')
ax31.set_ylabel('$dL/dt$')
ax31.semilogy(times, dldt)
ax31.tick_params(axis='both', which='both', labelsize=10)
ax31.grid(True)
ax31.semilogy(times, dldt, marker='o', ms=3, linestyle='None')
plt.savefig("%s_semilog_dldt.png" %params['timestamp'])
if __name__ == '__main__':
......
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